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Preface 


This book originates in a set of lectures I delivered at University Col- 
lege London between 2012 and 2018, corresponding to a master’s degree 
module in the Physics and Astronomy Department. Although a good num- 
ber of excellent published references on this material are available, my 
intended goal was to produce a set of notes that give a simplified, yet 
effective overview of the topic of galaxy formation and evolution, with 
special emphasis on dynamics. Extragalactic astrophysics started in earn- 
est as a discipline of physics when galaxies were discovered as “island 
universes” in the 1920s. The first steps in this field focused on under- 
standing the distance to the extragalactic nebula and to put them in 
context with the large-scale environment. Within the same decade, galax- 
ies became very distant objects and unveiled the cosmological process of 
expansion, shaping our current view of the Universe. Further analysis 
concentrated on the details of stellar dynamics in galaxies and the phys- 
ics driving the underlying components (dark matter, stellar populations, 
gas and dust). The advent of large galaxy surveys such as the Sloan Di- 
gital Sky Survey, and exquisite observations from facilities such as the 
Hubble Space Telescope have considerably transformed the field, allow- 
ing us to probe the distribution of galaxies over large cosmological scales, 
and to look at galaxy formation mechanisms ‘under the microscope’. The 
material presented in this book provides an introduction to the field for 
advanced undergraduates and beginning postgraduates. No substantial 
background in astrophysics is expected, but good knowledge of calculus 
is needed to enjoy the physics of galaxies at its fullest. 

I would like to thank the staff at UCL Press, especially Chris Penfold, 
for offering the opportunity to publish these lecture notes in an open ac- 
cess format. I would like to thank the previous lecturers of a precursor 
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module to the one I taught, Jonathan Rawlings, Jeremy Yates and Mark 
Cropper, for outlining and developing such an exciting course. I also 
thank Anna Pasquali, Prasenijt Saha, Witold Maciejewski, Andrew Hop- 
kins and Roger Davies for their support in putting together this book, and 
Joe Silk and Ofer Lahav for their guidance throughout the years. My inter- 
active audience, the students who participated in this module, are warmly 
thanked for their input, and often inquisitive minds. Iam especially grate- 
ful to Jennifer Chan, Lorne Whiteway and Ellis Owen. My wife, Isabel, will 
always have my immense thanks and gratitude for her undying support 
during the many weekends when these lecture notes were put together. 


Horsham, West Sussex, July 2018 
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1 
An introduction to galaxy formation 


Galaxies are the building blocks of the Cosmos. Separated by vast 
distances, they also serve as tracers of the cosmic expansion and the 
primordial density fluctuations that gave rise to structure in the universe. 
Galaxy formation requires an understanding of the most fundamental 
physical processes: gravitation, statistical mechanics, gas hydrodynam- 
ics, radiative transfer, atomic physics, etc. In this book we will focus on 
the gravitational side of galaxies, dealing with both the statistical treat- 
ment of galaxies as an N-body system evolving purely under gravitational 
forces and with the growth of galaxies from evolving density fluctuations 
in an expanding Universe. This introductory chapter presents an overview 
of the field, including the observables typically used to study galaxies, 
the mechanisms underpinning galaxy formation and the characteristic 
timescales involved. 


1.1 The main ingredients of a galaxy 


A galaxy is a complex system bound by gravity. In our current paradigm, 
the gravitational potential is dominated by dark matter, whose distribu- 
tion is much more extended than the visible part, and forms a spheroidal 
halo. The ordinary matter — loosely called “baryonic matter” — is made up 
mostly of hydrogen and helium, in the form of stars, diffuse and clumpy 
gas, dust, planets, etc. Although the dark matter dominates the mass 
budget - with a contribution of around 85 per cent in mass of the total 
matter content — emission in the electromagnetic spectrum is provided 
only by the baryons, except for potential, but hard to find dark matter 
particle annihilation events. Therefore, there is a substantial difference 
between mass and light in galaxies. 


The gaseous component provides fuel for star formation. A highly 
complex set of processes involving gas infall, turbulence, radiative trans- 
fer, feedback from star formation and magnetic fields plays a role in the 
physics of star formation (something we will leave aside in this textbook). 
In addition, dust provides an important tracer of star formation as it is 
typically found in gas-rich star-forming environments. The scattering, 
absorption and emission of radiation from dust makes this component 
key in the thermodynamics of star formation. Galaxies with very high 
star formation rates (starbursts) are often enshrouded in dust, with the 
most active regions being practically opaque to optical radiation and dis- 
playing prominent emission in the infrared by heated dust: this is the 
case with Ultra-Luminous Infrared Galaxies (ULIRGs) or submillimetre 
galaxies (SMGs). 

In addition to these components, it is worth noting the presence of a 
supermassive black hole (SMBH) at the centres of galaxies. With masses 
between a few million and several billion Suns, SMBHs can regulate the 
formation of their host galaxy. As gas accretes onto the SMBHs, a very 
luminous Active Galactic Nucleus (AGN) is formed. The energetic output 
from the AGN in the form of jets can affect star formation over the full 
scale of the galaxy, in ways that are still open to debate. 


1.2 Observables 


This section gives a nonexhaustive overview of the type of observables 
commonly applied to the study of galaxies. 


Colours 


In astrophysics, colour is defined as the flux ratio measured through dif- 
ferent filters (see section 1.5). The interpretation of a colour depends on 
the wavelengths covered by the filters. In the ultraviolet/optical/infrared 
spectral windows, colour can be considered a rough proxy of stellar age. 
Light from younger stellar populations is predominantly contributed by 
massive, luminous, blue stars. However, other factors — such as chemical 
composition or dust — will affect this interpretation: a red colour need not 
imply old stars. For instance, the red colours found in so-called ERO galax- 
ies (Extremely Red Objects) often originate from a young, but dusty stellar 
population. Figure 1.1 shows a mosaic of images of the nearby spiral 
galaxy M81, illustrating how a coverage of different regions of the elec- 
tromagnetic spectrum allows us to study different processes in galaxies. 
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X-ray © Witraviolet Optical lifrared . Radio (21cm) 


Figure 1.1 Different views of nearby galaxy M81 (NGC3031, distance 
3.7 Mpo). From left to right, images in the X-ray (NASA/CXC/Wisconsin/ 
Pooley & CfA/ Zezas), ultraviolet (NASA/JPL-Caltech/CfA/Huchra et al.), 
optical (NASA/ESA/CfA/Zezas), infrared (NASA/JPL-Caltech/CfA) and 
radio (NRAO/AUI/Adler & Westpfahl) spectral windows. 


The X-ray image reveals a diffuse component tracing hot gas; a central 
bright source betraying the presence of a supermassive black hole; and a 
number of point sources that correspond to X-ray binaries — stellar sys- 
tems where one of the members is a compact object (neutron star or 
black hole), whose strong gravitational potential drags and heats up the 
outer layers of the companion star. In contrast, the ultraviolet emission is 
due mainly to massive, young stars and reveals the sites of ongoing star 
formation. The optical and near infrared windows are dominated by the 
bulk of the stellar populations, whereas at longer wavelengths, in the far 
infrared, emission is produced by dust, that — like UV light — also traces 
sites of ongoing star formation. At even longer wavelengths, in the radio, 
emission originates from supernova remnants and HII regions (ionized 
hydrogen around star-forming sites), and at A = 21cm, we find resonant 
emission from neutral atomic hydrogen (HI). 


Spectroscopy 


The spectrum of a galaxy is the observed flux density as a function of 
wavelength, F(A), or frequency, F(v). Galaxy spectra carry valuable in- 
formation about the kinematics and the chemical composition of the 
stellar and gaseous components. Motions along the line of sight towards 
the observer affect the position and shape of the spectral features (both in 
absorption and emission) via the Doppler shift. For instance, the absorp- 
tion lines of a massive galaxy are significantly broader with respect to the 
same features in a low-mass galaxy, an effect caused by the higher velo- 
city dispersion of the stellar component. The bulk rotation of disc galaxies 
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Figure 1.2 Extragalactic spectroscopy at work: high quality stacked 
spectra from the Sloan Digital Sky Survey are shown for two early-type 
galaxies with different velocity dispersion — broadly tracing the mass of 
the galaxy. The insets zoom in special windows that feature absorption 
lines sensitive to the age, chemical composition and mass distribution 
of the underlying stellar component. (Source: data from Ferreras et al., 
2013, MNRAS, 429, L15.) 


is measured by the characteristic ‘S’-shaped pattern of the spectral line 
centres with respect to galactocentric distance. Moreover, a non-Gaussian 
analysis of the kinematic kernel — via higher order moments or a Gauss- 
Hermite expansion — allows us to further constrain the motion of the gas 
and stars of a galaxy. 

The absorption lines in ultraviolet/optical/infrared galaxy spectra 
originate in the atmospheres of their stars. Therefore, they carry inform- 
ation about the properties of the stellar populations (see an example in 
figure 1.2). For instance, spectral line strengths such as Balmer absorp- 
tion (Hf, Hy, H8), or metallicity-dependent features such as the Mgb-Fe 
complex in the 5100-5400A region provide constraints on the age and 
the chemical composition of galaxies. In emission, spectral lines originate 
from ionized regions, and the line luminosities constrain the proper- 
ties of the gas, including its composition, temperature and ionisation 
state. For instance, the BPT diagram! (named after Baldwin, Phillips and 
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Terlevich) provides a simple diagnostic to discriminate between emission 
from star-forming and AGN activity, by comparing the ratio between two 
pairs of line luminosities, such as H8/ [OIII] and Ha/[NII]. 

The dust component can also be probed with spectroscopy, includ- 
ing its overall attenuation effect with respect to wavelength and the pres- 
ence of spectral features, most notably the NUV 2175A “bump” or the 
silicate features in the infrared at 9.8 and 18um. Also in the infrared win- 
dow, emission lines from polyaromatic hydrocarbons (PAH) and thermal 
radiation are also sources of information regarding dust. 


Surface brightness 


Surface brightness (SB) is defined as the flux received from a section 
of the galaxy, i.e., within a solid angle. Therefore, it is not possible to 
measure the surface brightness of an unresolved source (e.g., the vast 
majority of stars observed with standard techniques). A simplified charac- 
terization involves the definition of elliptical isophotes, regions with the 
same surface brightness, leading to a one-dimensional surface brightness 
profile ©(@), where @ is the angular galactocentric distance, measured, 
e.g., as the semi-major axis of the ellipse describing the isophote. This 
distance can be translated into the physical projected two-dimensional, 
radius R, by use of the (angular diameter) distance (Da): R= Da tan 0 
(note R is often defined as a circularized radius vab, where a and b are 
the semi-major and semi-minor axis, respectively, of the corresponding 
ellipse that describes the galaxy). Although trivial for nearby sources, 
there is a second measure of the distance that translates between lumin- 
osity (L) and flux (F). This so-called luminosity distance (D;) is defined 
such that F =L/4rD?. For large but noncosmological separations (i.e., 
where the general relativistic effects are unimportant), Dg and D; are 
practically the same, making the surface brightness independent of dis- 
tance, but for cosmological separations D; = (1 + z)?Da, where z is the 
redshift to the source (see chapter 7). The radial distribution of the sur- 
face brightness of disc galaxies can be described by an exponentially 
decaying profile: 

ER) = Epe! Rh, (1.1) 


where h is the disc scale length, and £p is the SB at R=h. In contrast, 
elliptical galaxies have a much steeper profile (de Vaucouleurs, or R1/4 
profile): 


(R 1/4 
S (i) i (1.2) 
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The effective radius Re is defined such that half of the total flux is enclosed 
within Re, and Ee is the surface brightness at Re. A generic expres- 
sion often used in the description of the SB distribution is the Sérsic 
profile: 


1/n 
IR = Eee Le) - 1), (1.3) 


where n is the Sérsic index, and « = 1.9992n — 0.3271 is a normalization 
factor that ensures the flux within R < Rẹ is one half of the total flux.” This 
profile includes the exponential case forn = 1 and the de Vaucouleurs pro- 
file for n = 4. The value of the index is commonly used as a quantitative 
indicator of morphology, with the early-type galaxies having n > 2.5, and 
late-type galaxies having lower values of n. 


Morphology 


The morphology of a galaxy gives us information about the distribution of 
stars, gas and dust. Note that the appearance of a galaxy strongly depends 
on the spectral window (see figure 1.1). Morphological studies normally 
concern the stellar component as the main tracer of the gravitational po- 
tential of the galaxy. Galaxies are split roughly into three main groups: 
elliptical, spiral and irregular. The presence of a bar in spiral galaxies 
motivates a branching in the classification scheme, illustrated by the 
Hubble tuning fork diagram (figure 1.3). Elliptical and lenticular galaxies 
are combined into early-type galaxies and present a spheroidal distribu- 
tion, explained by a major merging event, or some collapse mechanism 
by which the total angular momentum was kept low. The oblateness of 
these systems cannot be fully explained by rotation (see chapter 5), and 
their spectral energy distribution corresponds to old and metal-rich stel- 
lar populations, which reflect an early, intense and efficient process of 
star formation and are corroborated by their chemical composition (see 
chapter 6). 

In contrast, spiral galaxies, more aptly described as disc (or late- 
type) galaxies, are flatter systems where a large fraction of the total kinetic 
energy is in the form of bulk rotation. For example, the (thin) disc of our 
Milky Way galaxy has a vertical extent about one-tenth of the disc size. 
The collapse of gas under gravity can develop such a rotating structure, 
as for instance, during the formation of our solar system. Disc galaxies 
have a more complex distribution of stellar populations than ellipticals, 
featuring ongoing star formation as well as a substantial presence of old 
stars. The central part of a disc galaxy usually hosts a spheroidal structure, 
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Figure 1.3 Hubble-de Vaucouleurs tuning fork diagram, showing the 
major morphological classification of galaxies into ellipticals (left) and 
spirals (right), the latter consisting of standard and barred spiral galax- 
ies. (Source: The images of the galaxies were created from observations 
taken by the Sloan Digital Sky Survey.) 


a bulge, whose origin also constitutes an open problem, with some bulges 
resembling an early-type galaxy (classic bulges) and others being the 
product of secular dynamical evolution (pseudo-bulges). Spiral arms are 
the most conspicuous features of disc galaxies. Their origin is based on 
dynamical resonances, which will be briefly explored in chapter 5. Irreg- 
ular galaxies are more complex dynamical structures, often betraying the 
presence of an ongoing merger or tidal interaction with a neighbour. 

The standard method of morphological analysis involves visual in- 
spection of images, preferably done through several filters. Morphology 
can also be determined from the surface brightness profile (see above), 
or by the application of alternative methods involving nonparametric ob- 
servables, which do not make any assumption about the radial surface 
brightness profile. A number of observables are defined such as the con- 
centration, asymmetry or clumpiness, or even higher order moments of 
the pixellated surface brightness distribution. More recently, machine 
learning methods are being applied to perform “visual classification” in 
a fully automated way, using, for instance, artificial neural networks 
trained on visually classified data sets.* 
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Size 

It is difficult to provide a clear-cut definition of the size of galaxies. 
Being diffuse objects, their borders are fuzzy. One typically simplifies the 
problem to a one-dimensional equivalent by binning the observation radi- 
ally within elliptical isophotes. A traditional definition of size, D25, hypo- 
thesizes that the galaxy extends in a region brighter than 25 mag arcsec” 
in the B band. This choice is motivated by the fact that surface bright- 
ness does not vary with distance over noncosmological scales. However, 
we know that over large distances, the surface brightness does indeed 
change, decreasing as (1+2)*, where z is the redshift (termed ‘Tolman 
dimming’). A more robust criterion is based on the Petrosian radius, 
derived from the following expression: 


=x (Ro) 


ee 1.4 
(Z (R))R<Ro ee 


n(Ro) = 


This function starts at n ~ 1 for Ro ~ O and decreases outwards. The Pet- 
rosian radius is defined as the value Rp for which n(Rp) = 0.2. Its being a 
ratio of surface brightness eliminates the dependence on the cosmological 
Tolman dimming. Another option involves fitting the observed data with 
a generic surface brightness profile such as the Sérsic function presen- 
ted above (equation 1.3), so that the parametric effective radius of the 
function can be used as a measure of size. 


Exercise 1.1 


Find the Petrosian radius as a function of h for the exponentially 
decaying surface brightness case shown in equation 1.1. 


Luminosity function 


Galaxies can be classified according to their absolute luminosity (i.e., 
their power, or energy emitted per unit time, usually defined with respect 
to solar luminosity, Lo, or given by their apparent flux when located at a 
fiducial distance). The volume number density of galaxies (n) per lumin- 
osity interval is the Luminosity Function. It can be suitably described by 
a power law with an exponential cutoff at the bright end, defined as the 
Schechter function: 
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onal = (L)dL= oo( g) e =, (1.5) 
with three free parameters: the characteristic luminosity (L,, usually 
given in units of Lo, or as an absolute magnitude, M,), the power law 
index, equivalent to the slope at the faint end (œ), and the normalisation, 
given by the number density of L, galaxies (given, aside from a factor e, by 
®o, usually in units of Mpc~?). We will see in chapter 7 that the shape of 
the Schechter function is a direct consequence of galaxies being created 
from random density fluctuations following collapse once the fluctuations 
are large enough. 

The parameters describing the luminosity function change with 
respect to galaxy morphology. For instance, elliptical galaxies are less 
numerous (lower ®o), but are brighter (higher L») than disc galaxies. 
Moreover, these populations change with respect to environment, with 
ellipticals being more prevalent in galaxy clusters (a property often 
termed the ‘morphology-density relation’). An alternative description of 
the galaxy census is the stellar mass function, with a similar notation as 
in the Schechter function, but referring to mass instead of luminosity. We 
will see below that a comparison between the stellar mass function and 
the dark matter mass function expected from a simple theoretical argu- 
ment gives a very powerful diagnostic of the nontrivial relation between 
dark matter halos and galaxies. 


Exercise 1.2 


Show that if we adopt the Schechter luminosity function (equa- 
tion 1.5), we need a >-—2 if the total luminosity (integrated 
throughout the whole galaxy population) is finite. Then consider 
the œ =—1 case: find the total luminosity from all galaxies, and 
show that the total number of galaxies diverges. How can an infinite 
number of galaxies have a finite total luminosity? 


Star formation rate 


The star formation rate (SFR, often denoted by the greek letter y) is 
defined as the mass content in stars created per unit time. This can be 
defined over different scales as an integrated quantity over a cosmological 
volume, for a single galaxy, or as a volume (or projected surface) density 
of star formation within a galaxy. The average star formation rate of 
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the Milky Way galaxy at present is ymw ~ 1Mo yr~!, and star-bursting 
galaxies can sustain rates hundreds (or even thousands) of times higher. 
Measurements of the SFR use indirect tracers, focusing either on the 
presence of young stars or on the fuel that triggers formation. The ultra- 
violet emission from a galaxy reveals the presence of massive (therefore 
young) stars. Emission lines from ionized gas also betray the presence 
of the most massive stars, hot enough to ionize large volumes surround- 
ing the star-forming regions. The (hydrogen) gas can be observed both 
in atomic and molecular form (the latter normally detected indirectly 
through the presence of CO), and the dust — which is abundant in star- 
forming regions — can also be probed by attenuation of starlight or by 
emission at infrared wavelengths. There is a whole industry of extracting 
star formation rates from observations, beyond the scope of these lecture 
notes.” 


Environment 


‘Environment’ refers to the large-scale regions where galaxies are found. 
It can be characterized by a density averaged over volumes that include 
many galaxies. It provides a link between local formation processes 
in galaxies and the large-scale structure of the Universe (figure 1.4). 
Roughly four different environments can be considered: voids (under- 
dense regions, away from large-scale filaments), field (representing aver- 
age densities over cosmological volumes), groups (gravitationally bound 
structures that include several massive galaxies, with velocity dispersions 
among group members comparable with the velocity dispersion of stars 
within galaxies, i.e., 100-300 kms~') and clusters (the highest density 
regions, including thousands of galaxies within a ~1 Mpc radius, and ve- 
locity dispersions >500 kms~!; see chapter 9). Another classification of 
environment can be given by the mass of the dark matter halo that hosts 
a given galaxy. Moreover, inside groups and clusters, the environment 
of the central region is very different from the cluster outskirts, and the 
distinction between a central (i.e., the most massive galaxy in a dark mat- 
ter halo) and a satellite appears to be fundamental in the description of 
galaxy formation.° 

There are important environment-related processes that affect the 
evolution of galaxies: mergers, strangulation, ram pressure stripping, har- 
assment. The environment is tightly linked with the initial conditions of 
the system: clusters are structures formed from higher (and rarer) fluctu- 
ations of the primordial density distribution. Furthermore, because of the 
inherent growth of structures with cosmic time (bottom-up formation), 
clusters can be found only at relatively later cosmic epochs. 
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Figure 1.4 A slice of the Universe as seen by the EAGLE galaxy form- 
ation project. The cosmic web is readily apparent, with filaments, 

voids and clusters (at the confluence of filaments). The zoomed-in 
panels give an idea of the size of a typical galaxy within the cosmic 

web. Environment-related processes will shape the formation histor- 

ies of galaxies. (Source: Schaye et al., 2015, MNRAS, 446, 521.) Used by 
permission of Oxford University Press. 


Nuclear activity 


Most galaxies harbour a supermassive black hole (SMBH) at their centres, 
with a mass in the range 10°-10?Mo. When gas is available, as it flows 
into the SMBH, it forms an accretion disc where the material gets heated, 
emitting copious amounts of radiation across a wide range of wavelengths 
and producing an AGN that can outshine the whole galaxy. Only a few 
galaxies display AGN activity, but it is believed that most massive galaxies 
experienced AGN activity some time during their past formation history. 
There are many types of AGN: Quasars, Seyferts, Radio Galaxies, Liners, 
BL Lacs, etc. 

The fundamental issue related to galaxy formation is that AGN 
activity can affect the flow and heating of the gas that fuels star formation 
in galaxies. Therefore, although AGNs are confined to a minuscule volume 
within the galaxy (a 10?Mo SMBH has a radius rg = 2GM/c* ~ 20 AU 
and a region of influence r; = GM, /o* ~ 50 pc), they can control galaxy 
formation (i.e., over scales 250 kpc). It is believed that AGN activity 
provides a feedback mechanism that contributes to the quenching of star 
formation in massive galaxies. A strong correlation is found between 
galaxy properties and the mass of the central black hole.” 
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Distance 


A simple but fundamental observable is the distance to galaxies from 
our vantage point. Measuring distances beyond the solar system entails 
a progressive set of steps (called the ‘cosmological distance ladder’) to 
measure successively larger distances. The first ‘extragalactic’ rung of the 
ladder is based on the use of variable stars, especially bright Cepheids. 
Their period-luminosity relation allows us to determine their true lumin- 
osity (L) from a sequence of observations with time (the light curve). In 
combination with the observed flux (®), the distance (d) is obtained via 
® = L/4rd?. This method pioneered the field of extragalactic astrophys- 
ics, when our neighbour, the Andromeda galaxy, was found to lie much 
farther out than any of the visible stars. Additional steps of the ladder 
involve the use of tight scaling relations among galaxies, such as the Tully- 
Fisher relation (chapter 5). Closer to home, the first step of the ladder 
uses basic trigonometry: by carefully measuring the positions of stars, 
we find that in addition to their proper motion (caused by their relat- 
ive velocity with respect to the Sun), they follow a conspicuous yearly 
cycle, with all stars moving “in phase” but with different amplitudes. 
This motion is caused by parallax, the apparent change in the observed 
position as the Earth orbits the Sun. Therefore, those with the largest 
cycles are closest to us. These cycles are nevertheless very small, and 
thus limit useful estimates of parallax to the nearest stars. A natural 
unit of distance is therefore the one at which a star would have a paral- 
lax of 1 arcsec. This is called the ‘parsec’ (pc), and trivially amounts to 
~1AU/z, where z is the parallax in radians and AU denotes an astro- 
nomical unit (the average Sun-Earth separation, about 150 million kilo- 
metres). Therefore, 1 pc amounts to 3.086 x 101° m, or 3.26 light-years. 
Distances within galaxies are measured in kiloparsec, and galaxy-galaxy 
separations are measured in Megaparsec (Mpc) or Gigaparsec (Gpc). The 
distance from the Sun to the centre of the Milky Way galaxy is approx- 
imately 8 kpc, and the distance to the neighbouring Andromeda galaxy 
is 0.8 Mpc. 


Redshift 


At the top of the cosmological distance ladder, the separation between 
galaxies takes advantage of redshift. Due to the expansion of the Universe, 
photons with wavelength ào, emitted from a distant galaxy (at cosmic 
time t), are observed (at time tọ), at a longer wavelength, obs, stretched 
by a factor z, termed ‘redshift’: 
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À t 
obs ga 0) 


ho ~ a(t)’ 


(1.6) 


where zis the redshift, and a(t) is the scale factor of the Universe, a dimen- 
sionless quantity that represents the global cosmic expansion. The scale 
factor at present time, a(to) can be arbitrarily set to one. Since redshift 
increases towards earlier times, it directly correlates with distance. Note 
that in groups of nearby galaxies the relative motions (termed ‘peculiar 
velocities’) will also affect, via the Doppler effect, the observed redshift. 
Only at very large distances does the cosmological effect overwhelm any 
contribution from peculiar velocities. We will see in chapter 7 that it is 
convenient to separate the uniform expansion (termed the ‘Hubble flow’) 
and the peculiar motion, by use of co-moving distances between galax- 
ies where the contribution from the uniform expansion is factored out. 
We will also see in that chapter that the concept of distance becomes less 
trivial, and different distances to the same object will be obtained when 
using alternative definitions (section 7.1). Measuring galaxy properties 
over a range of redshifts implies testing the evolution of galaxies across 
cosmic time. As reference, at redshift z=2, galaxies are found in a Uni- 
verse that is only 3 Gyr old, i.e., about 20 per cent of its current age. At very 
high redshift, cosmic age poses a stringent constraint regarding the times- 
cales expected in the formation of galaxies. For instance, bright quasars 
or galaxies at z~10 have only about 500 Myr of time to evolve, so obser- 
vations at such high redshift substantially reduce the number of potential 
scenarios of formation. 


1.3 Physical processes 


We switch from the observational properties of galaxies to the physical 
mechanisms of galaxy formation. This section gives a succinct view of 
some of the processes that play an important role in the formation and 
evolution of galaxies. 


Initial conditions 


Galaxies develop from small density perturbations in the dark matter 
distribution at early times. Under gravity, these perturbations grow 
within an expanding Universe, until virialisation is achieved, forming 
stable structures called ‘dark matter halos’. In the current cosmological 
framework, about one-sixth of all matter is in the form of baryons (mostly 
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H and He). Originally in gaseous form, this material dissipates energy — 
in contrast with the collisionless dark matter. The gas falls to the centre of 
the halo, cooling down and eventually forming stars. 

Certain initial conditions in the overdensities are needed to ex- 
plain the properties of galaxies: too low and there will not be enough 
galaxies at present; too high, and we would be surrounded bya large num- 
ber of massive, old galaxies, at odds with the observational evidence. 
The present distribution of galaxies is compatible with a model where 
~24 per cent of the total matter/energy content of the Universe is in 
the form of dark matter, an additional ~4 per cent is in the form of ba- 
ryons and the remaining ~72 per cent is made up of an unknown energy 
field commonly known as ‘dark energy’, that produces an accelerated ex- 
pansion (this is the standard ACDM framework, where A refers to the 
dark energy term, and CDM stands for cold dark matter). Galaxy form- 
ation is inherently tied to cosmology, and moreover to the fundamental 
properties of particles. In fact, many cosmology-oriented experiments 
use large galaxy surveys to map the large-scale structure of the Universe 
(e.g., the Sloan Digital Sky Survey, the Dark Energy Survey, or the Euclid 
mission). 


Gravitational instability 


Small fluctuations are rapidly amplified by gravity — as neither expansion 
nor the repulsive effect of dark energy are effective over galaxy scales. 
It is useful to quantify the evolution of the fluctuations with the density 
contrast, defined as follows: 


PĒO 
(ot) 


SË: = (1.7) 


In a noncosmological environment, where the effects of expansion are 
negligible (e.g., a gas cloud within the Milky Way), a small overdensity 
will grow exponentially with time (5 œ et). Over cosmological scales, the 
expansion will slow down this growth rate to a power law of cosmic time, 
ôx t”, with the power law index depending precisely on the expansion 
rate, given by the cosmological model. Cosmologists often work in the 
linear regime (6 < 1), where the growth equations can be expanded as a 
Taylor series, and truncated to the lowest order (i.e., linear) terms. In con- 
trast, astrophysicists tend to work instead in the highly nonlinear regime 
(6 >> 1), where galaxies constitute huge overdensities with respect to the 
background density of the Universe. 
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The hierarchical growth of structure inherent to the current cos- 
mological paradigm implies that fluctuations over small scales are the 
first to form virialized, stable systems at early times, whereas the more 
massive structures assemble later. Therefore, the linear regime corres- 
ponds at present time to very large scales (of cluster/supercluster size), 
and the equivalent linear regime for galaxy scales took place at earlier 
times. Hence, cosmology studies using galaxies need to probe very 
large volumes. Over smaller scales (galaxy clusters or groups), galaxy- 
related processes overwhelm any imprint on the observations related to 
cosmology. 


Gas cooling and star formation 


The only channel of interaction among dark-matter particles is via 
gravity, in a collisionless way, i.e., simply following the trajectories dic- 
tated by the overall gravitational potential. In contrast, gas particles 
(hydrogen/helium) can lose energy, as they also interact through elec- 
tromagnetic forces. Photons emitted from collisions of gas particles 
are lost from the system, and act as an energy sink. The cooling rate 
depends on the physical properties of the gas: temperature, density, 
composition. At very high temperatures (T > 107K), the gas is fully 
ionized, and cools mainly through bremsstrahlung emission from free 
electrons. At lower temperatures (104 K < T <10” K) atomic recombina- 
tion/collisional excitation processes dominate, in ways that are strongly 
dependent on the species present (i.e., the metallicity of the gas). At 
even lower temperatures (T < 104 K), the gas is almost completely neut- 
ral, and collisional excitation of fine structure energy levels (atoms/ions) 
and rotational/vibrational energy levels (molecules) will contribute sig- 
nificantly in this regime. We will see in chapter 7 that the balance between 
gas cooling and the dynamical timescale from gravity is responsible for 
the characteristic sizes and masses of galaxies. To quantify the role of 
cooling, it is convenient to define a cooling timescale as the time it would 
take to remove the internal energy of the gas at the rate dictated by its 
cooling. It is proportional to (T/nA), where T is the temperature of the 
gas, n is its number density, and A is the cooling rate, a complex func- 
tion that, in its simplest form, depends on temperature and chemical 
composition. 

The cooling of atomic gas leads to pressures and temperatures at 
which molecular gas can form. Following the formation of molecular 
clouds, fragmentation of clumps create the so-called pre-stellar cores, 
which act as seeds of star formation. This process poses one of the 
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most difficult and challenging problems in physics. For instance, the 
distribution of stellar masses in a newly forming region (called the ‘Ini- 
tial Mass Function’, see chapter 6) remains an open problem as analytic 
models cannot describe the complexity of the underlying processes, and 
hydrodynamical simulations lack the necessary dynamical range to probe 
these processes in detail. 

From a phenomenological point of view, we make use of scaling rela- 
tions backed by the observations that allow us to describe in simple terms 
the process of star formation in galaxies. One of the most fundamental 
scaling relations is the one between the star formation rate (SFR) and the 
gas content. The SFR is found to follow a power law of either the volume 
gas density (Schmidt law) or the surface gas density (Kennicutt law). For 
instance, the latter is defined as follows: 


= 1.6 x 10-#(— =e r (1.8) 
1Mọopc-2 

The power law index can be derived from a straightforward theoretical 

argument: Y & pgas/tdyn, Where tdyn X paas" 2 is the dynamical timescale 

of collapse of the gas (see exercise 1.3). More recent observations sug- 

gest this law is valid when the densities correspond to the molecular (not 

atomic) gas.” 

Figure 1.5 illustrates the overall link between the mass of a galaxy 
and its formation history. The left panel shows the colour-mass diagram 
of low redshift galaxies — colour can be used as a rough proxy for age, 
with red colours (g—r= 0.8) meaning old populations. At the massive 
end, most galaxies — at late cosmological times — populate an elongated 
region towards the red side (top of the figure) called the ‘red sequence’, 
consisting of old stars. These are galaxies whose stellar component was 
formed at early times. In contrast, low-mass galaxies are markedly blue, 
populating a region called the ‘blue cloud’ and representing star-forming 
systems. The panels on the right give a very schematic view of the inter- 
play between the red sequence and the blue cloud. Mergers will induce 
mass growth in galaxies. The ages of the stellar component reveal the age 
of the progenitors that merge into a more massive system, as well as poten- 
tial star formation, when the merger includes gas (black arrows, termed 
‘wet mergers’). Mergers between galaxies that have no gas (white arrows, 
‘dry mergers’) mix their stellar content without adding any younger com- 
ponent. However, this picture is, alas, too simplistic. Incidentally, the 
region between the blue cloud and the red sequence is called the ‘green 
valley’, and it may hold valuable clues about galaxy evolution processes, 
as it represents a transition stage from star-forming galaxies to quiescent 
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Figure 1.5 Evolution of galaxies on the colour versus mass diagram. 
The left panel shows the observed distribution of galaxies from the Sloan 
Digital Sky Survey (https://www.sdss.org), coding the number dens- 
ity of galaxies as a grey scale. The panels on the right give a number of 
schematic trends of evolution, with black (white) arrows representing 
gas-rich (gas-poor) mergers. (Source: adapted from Faber et al., 2007, 
ApJ, 665, 265.) 


systems. It is in the green valley where we expect the feedback-driven 
processes that lead quenching to operate. 


Feedback 


The simple laws that relate the star formation activity in a galaxy to 
the amount of gas — such as the Kennicutt law of equation 1.8 — also 
depend on additional properties that can be encoded as a star forma- 
tion efficiency. For instance, in its basic form, Schmidt’s law (y = voga) 
leaves a free parameter (v) that relates to the efficiency of the process. 
Various mechanisms can reduce the star formation rate, even if gas is 
present. The most direct one is star formation itself: massive stars end 
their lives as supernova explosions, releasing vast amounts of energy both 
as a heat source and as ejected material ploughing through the gaseous 
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Cosmological Limit 
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Figure 1.6 The ratio between galaxy (stellar) mass and dark matter 
halo mass is shown as a function of galaxy mass, at two different red- 
shifts. The top horizontal line is the cosmological limit (~ i) derived 
from the Planck 2015 cosmological parameters (see table 7.1). 
(Source: adapted from Behroozi et al., 2010, ApJ, 717, 379.) 


component that feeds star formation. Since star formation is a highly 
clustered process, many supernova explosions can trigger a super-wind 
that will drag substantial fractions of the gas into intergalactic space, in 
a form of negative feedback (termed ‘stellar feedback’). A second feed- 
back process involves the activity from a central Active Galactic Nucleus, 
ejecting gas and heating the surrounding regions, in a similar way to the 
above, and suppressing star formation. 

One of the key figures that conveys the most fundamental issues of 
galaxy formation (figure 1.6) compares the mass function of dark matter 
halos (from numerical simulations of structure formation) with the ob- 
served stellar mass function of galaxies. Making the assumption that one 
can match the abundance of dark matter halos to galaxies (the former be- 
ing hosts to the latter), it is possible to define the mass ratio between dark 
matter halos and stellar mass in galaxies (this procedure is called ‘abund- 
ance matching’). The figure shows that at both ends (massive galaxies 
on the right and dwarf galaxies on the left), the mass ratio decreases, 
suggesting that a mechanism must be present to reduce the efficiency 
of conversion from gas into stars. Stellar feedback is believed to con- 
trol the low-mass end, producing an increasing correlation of the ratio 
with galaxy mass: lower mass galaxies will have shallower gravitational 
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potentials, ejecting more gas, and therefore reducing the star formation 
rate. In contrast, massive galaxies have stronger potential wells, mak- 
ing stellar feedback less effective. AGN-driven feedback is supposed to 
drive the negative correlation at the high-mass end, as the mass of the 
SMBH - and therefore the energetic output from the AGN — will increase 
with galaxy mass. Also note that at the peak of the trend, correspond- 
ing to galaxies with stellar masses ~3x 10!°Mpo, the ratio is about 4 per 
cent, substantially lower than the ‘primordial’ ratio between baryons and 
dark matter (~16 per cent, top horizontal line). Moreover, galaxies at the 
peak of the ratio do not have a vast amount of gas (or dust) to compensate 
for this mismatch, suggesting that only a small fraction of the available 
baryons is found in galaxies, with the remainder either not forming part 
of the galaxy formation process or being ejected during the main forma- 
tion phases. The redshift evolution, from the dashed line (at z = 1) to the 
solid line (at z= 0.1) illustrates the significant growth in the dark mat- 
ter halo at the massive end, without a comparative stellar mass growth — 
another sign of the mismatch between the growth of dark matter halos 
and galaxies. 


Galaxy mergers and dynamical evolution 


Dark matter halos (and their galaxies) interact among themselves, lead- 
ing to a hierarchical growth process, whereby small structures — formed at 
early times — lead to more massive structures via mergers. When halos of 
a similar mass merge (major mergers), the gas in the progenitors is shock- 
heated to the new equilibrium (virial) temperature of the newly formed 
halo. Galaxies in the progenitor halos will also merge, producing a new 
galaxy at the centre of the newly formed halo. When the mass ratio of 
the halos is larger (minor merging), the merging may not be so effective, 
keeping the galaxy from the small halo as a satellite. Merging processes 
lead to morphological changes that strongly depend on the environment. 
For instance, spheroidal galaxies, believed to form from major merging 
processes, are more abundant in the higher density regions of clusters 
than in the field. Massive galaxies are more compact at high redshift, 
suggesting a substantial growth in size through galaxy mergers. 

Internal dynamical effects may also play a role in the evolution 
of a galaxy. Instabilities can redistribute mass and angular momentum, 
changing the morphology of a galaxy. For instance, a thin disc with 
a high surface mass density is susceptible to instabilities that produce 
structures similar to those found in barred spiral galaxies. This process 
will efficiently drive gas towards the centre of the galaxy, fuelling star 
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formation or AGN activity. Once formed, a bar can buckle away from the 
plane of the disc, to form a spheroidal object, called a ‘pseudo-bulge’,!° 
which should be distinguished from classic bulges, themselves formed 
via mergers. Hence, the morphology of a galaxy can be driven either by 
external mechanisms such as the merger history or by internal (secular) 
processes. 

Dynamical evolution introduces a characteristic timescale. In 
exercise 1.3 we find that the collapse time of a gravitating system with 
negligible internal kinetic energy is uniquely determined by its dens- 
ity (assuming no other forces oppose gravity). The dynamical timescale 
(tayn X p712) broadly introduces a lower limit on the duration of the 
‘settling’ time that a galaxy undergoes after a dynamical event such as a 
merger or sudden mass loss. 


Exercise 1.3 


The dynamical timescale can be defined by the time it would take 
for a system driven only by gravity to collapse to a point (neg- 
lecting the initial contribution from kinetic energy). We assume 
each particle moves only under the gravitational potential of the 
collapsing mass, neglecting any nongravitational effects or local in- 
teractions with neighbouring particles. Show that this timescale 
depends only on the density of the galaxy as follows: 


TRE 370 
dm 41 32Gp' 


Note that for the Sun (pọ = 1.4 g cm™°?), this time is approximately 
30 min, i.e., the time it would take to collapse if it suddenly trans- 
formed into pressureless dust. This type of problem will reappear 
in chapter 7 with Friedmann’s equations and the spherical collapse 
model. 


Chemical enrichment 


The ‘primordial’ periodic table of the elements is comprised only of hy- 
drogen (76 per cent), helium (24 per cent) and traces of lithium. The 
phase of cosmological nucleosynthesis — when the temperature and dens- 
ity are appropriate for an efficient processing of nuclei — happened a few 
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minutes after the Big Bang, but it was not long-lived enough to allow for 
the efficient creation of nuclei more massive than helium. The chemical 
composition of galaxies is determined, instead, by stellar nucleosynthesis. 
As stars evolve, new elements are synthesized in their hot and dense in- 
teriors via either thermonuclear reactions in the core, neutron capture, or 
explosive processes in the most massive stars. These elements are released 
during the final stages of stellar evolution, contributing to the change in 
composition (metallicity) of new generations of stars. 

Hence, the evolution of the composition of stars with age (chem- 
ical enrichment) constitutes an important tracer of the formation history 
of a galaxy, as we will see in chapter 6. For instance, we can use the 
metallicity of stars as a rough proxy of age: extremely metal-poor stars 
are relics from the early formation phases of the Milky Way. However, the 
efficiency of star formation can alter this result: the cores of massive ellipt- 
ical galaxies feature high metallicity, but they are also very old, reflecting 
a very intense process of formation that quickly enriched the interstel- 
lar medium at early times. The positive correlation between galaxy mass 
and its metallicity (measured both in the stellar populations and in the 
gaseous component) can be explained via feedback processes that trigger 
a metal-rich wind, preferentially removing material from the shallower 
gravitational potential wells of lower mass galaxies. 


Evolution of stellar populations 


The stellar populations in galaxies feature a wide range of stellar mass, 
age, and chemical composition. Observational data in the ultraviolet/ 
optical/infrared spectral window, covering approximately a wavelength 
interval from 0.1 to 4 micron, allow us to study the underlying stellar pop- 
ulations. Synthetic models combine our knowledge about stellar structure 
and evolution, producing spectra for a choice of parameters (age, metalli- 
city, initial stellar mass distribution, etc.).!! These models can be used to 
backtrack the star formation histories of galaxies. Although fraught with 
degeneracies inherent to the complexity of the population mixtures and 
the severe overlapping of spectral features on the stellar atmospheres, 
these models have allowed us to shed light on the connection between 
various observational properties and the star formation histories of galax- 
ies. For instance, it has been found that stellar mass is the dominant 
driver of the overall formation history of elliptical galaxies, with massive 
galaxies featuring an earlier and more intense process of formation, 
whereas low-mass galaxies extend their formation over much longer 
timescales. 1? 
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1.4 Stellar clusters 


Smaller than galaxies, stellar clusters are associations of stars with total 
stellar mass below ~10°M.. They constitute sites of recent or past star 
formation — watch out for a notation overload, as we also use the term 
‘clusters’ to refer to galaxy clusters! The population of stellar clusters is 
split into two main categories. (1) Globular clusters (GCs) are highly con- 
centrated, spherical stellar systems whose trajectories do not follow the 
disc. Dominated by old stellar populations with relatively low metallicit- 
ies, GCs are formed during an early and intense period of star formation. 
(2) Open clusters (OCs) are stellar associations found in the disc of the 
galaxy. They represent sites of recent or ongoing star formation, and their 
chemical composition is more evolved than that of GCs. As a stellar cluster 
ages, the system dissolves into the general distribution of disc stars via 
various mechanisms, induced mainly by tidal forces within the galaxy. 
OCs are less massive than globular clusters (with stellar masses around a 
few x10°Mo). GCs and OCs in our galaxy and nearby systems (most not- 
ably the Magellanic Clouds) can be observed in detail. In distant galaxies 
it is possible to observe GCs as unresolved sources, and they can be used 
as dynamical tracers of the underlying gravitational potential. Moreover, 
star-bursting galaxies feature super-star clusters, very massive and young 
stellar associations that can be considered a young version of GCs. 


1.5 A technical note on astronomical observations 


Flux density (or irradiance) is the main observable of an astronomical 
source. It is defined as the energy received per unit time and surface (as 
in, e.g., the collection area of a telescope). Typical units are erg s~! cm~? 
(cgs). However, flux densities are often quoted as magnitudes on a logar- 
ithmic scale: 


m=Z -— 2.5logF =Z — 2.5 log (1.9) 


L 
4nD2’ 
where Z = 2.5 log Fo is the zero point, which corresponds to the flux of 
a reference star that, by definition, has zero magnitude. Traditionally, 
Vega (« Lyr), one of the bright stars in the (Northern) Summer triangle 
has been used as reference (see below). D is the distance to the source, 
and L is the luminosity (or power, energy emitted per unit time). In 
the above equation, m is the apparent magnitude. In order to compare 
the luminosities from sources located at difference distances, it is use- 
ful to factor out the distance by defining the absolute magnitude (M) as 
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the apparent magnitude that the same source would have if located at 
a fixed, reference distance (dg =10 pc). Since the flux decays as 1 /D?, 
we have 
ere (1.10) 
do 


The term 5logD/do is called the ‘Distance Modulus’. This is a simple 
expression where redshift and other cosmology-related effects are neg- 
ligible. So far we did not take into account the fact that a source 
emits photons over a range of wavelengths, according to some dis- 
tribution (called the ‘spectrum’). We define F(A) as the observed flux 
density per unit wavelength, also called the ‘spectral irradiance’ (meas- 
ured in, e.g., ergs-!cm~2A~!). Alternatively, the spectrum can be 
given per unit frequency, F(v); a standard unit is the Jansky (1 Jansky= 
10-3 ergs~! cmt Hz~!). Note that F(v) Æ F(A), as one is measured per 
unit frequency and the other is measured per unit wavelength. We can 
relate both using dF = F(A)di = F(v)dv, so that 


2 
F(v) = FQ). (1.11) 


The total, or bolometric, flux consists of energy from all wavelengths, or 
frequencies: F = f F(A)da = f F(v)dv. In spectroscopy, the observations 
are given as either F(A) or F(v). Another expression often used to describe 
the spectrum is AF(A), proportional to the number density of photons 
received on a detector per unit time. In photometry, we usually define 
the flux density passing through a filter with a transmission throughput 
S(A) as 
= So" FO)SQ)dX Ea À (1.12) 

Jo SOA)dA 
from which we can determine the filter-specific apparent magnitudes, as 
follows: 

ms = Zs — 2.5 log Fs. (1.13) 


For example, in the ultraviolet/optical/infrared windows, there are sets 
of standard filters, such as those shown in table 1.1. The choice of ref- 
erence star — which gives the ‘zero points’ Zs — defines the photometric 
system used. In each one, the reference standard has, by definition, zero 
magnitude through any passband. The most standard systems are: 


1. VEGAMAG, where the reference star is Vega (aLyr); 


2. AB, where the reference star has F(v) = 3, 631 Jansky, constant at all 
frequencies; 
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Table 1.1 Some of the standard photometric passbands, showing the central 
wavelength (Ao) and the full width at half maximum (AA) of the transmission 
function. For reference, we include the absolute magnitude of the Sun (Mo). 
Note that the apparent magnitude of the Sun from the Earth’s position in the 
V band is Vo = — 26.74. 


Filter U B V R I J H K 


Ao (A) 3650 4400 5500 7000 9000 12500 16500 22000 
AA CÅ) 680 980 890 2200 2400 3800 4800 7000 
Mo +5.61 +5.48 44.83 +4.42 44.08 +3.64 +3.32 +3.28 


3. STMAG, for which the reference star has F(A) = 3.631 x 107? ergs"! 
cm~? Å-1, constant at all wavelengths. 

A colour index (or ‘colour’, in short) is defined as the ratio of flux densities 
of the same source, when observed through two different filters. In the 
logarithmic magnitude scale, the colour is therefore the difference of the 
magnitude in each of the filters; for instance, 


0 


L L 
U—V=my -my =My My =2.5(log 7< -log =r), (1.14) 
U 


where the 0 superindex refers to the reference star in a given photomet- 
ric system. We use luminosities in this expression because the ratios of 
fluxes and luminosities should be the same. The colours of astrophysical 
sources provide useful information about their properties: the temperat- 
ure of a star, the age and chemical composition of a stellar population, or 
the redshift of a distant galaxy. Colours contrasting different wavelength 
regions are sensitive to specific properties of the stellar populations or 
the gas and dust content. We can interpret a set of photometric colours 
as the low-resolution equivalent of a spectrum. Indeed, many studies 
are based on fluxes covering a large set of passbands that jointly provide 
a (multicolour) spectral energy distribution. The spectral resolution is 
R=A/Ad, where å is the (average) wavelength of the observation and 
Ad is the width of the passband (commonly quoted as the full width 
at half maximum of the spectral response function). From table 1.1 we 
find that typical broadband photometric studies have an effective resol- 
ution around R~ 5. This is equivalent to smoothing the spectra shown 
in figure 1.2, reducing each spectrum to just five fluxes within the same 
wavelength interval. Although spectroscopy provides much more detailed 
information about galaxies, broadband photometry allows us to probe 
deeper, as we collect all the photons within the response of the filter into 
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a single ‘bin’. The standard approach in observational work is to survey 
and select interesting targets with photometry, following up a targeted 
set with spectroscopy, where the typical spectral resolution adopted when 
observing galaxies lies between R = 500 and 5,000. As acompromise, one 
can use medium-band or even narrow-band filters to increase the spec- 
tral resolution of photometric measurements (no higher than about 100 
or 200) but allowing us to reach deeper than with spectroscopic obser- 
vations. Note that the dynamical state of galaxies introduces a velocity 
dispersion of the stellar component o > 50kms~!. Therefore, a spectral 
resolution above R=c/o ~ 6, 000 is less useful for the analysis of, say, ab- 
sorption line features originating in stellar atmospheres. Such a natural 
resolution limit imposed by the motion of stars in galaxies complicates, 
for instance, the study of detailed line strengths, which appear heavily 
blended in the spectra of galaxies with respect to data from individual 
stars. 


Exercise 1.4 


Consider two filters with a top-hat spectral response as shown in 
the figure. Compute the colour S1 — S2 for a star with a spectrum 
F(A) =k, where k is a constant. Use the AB photometric system. 


S,(A) 


S,(A) 


m. 


Throughput 
o 
on 


[o) 


1 2 3 1 2 3 
Wavelength (A/A,) 
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z 
The classical theory of gravitation 


The Newtonian theory of gravitation is presented here, along with the 
most fundamental equations. We will begin with Kepler’s problem, the 
simplest case of a gravitational system. The centrepiece of classical gravity 
as a field is Poissons equation, which relates the mass density distri- 
bution with a scalar function (potential) from which the forces (and 
hence the dynamical evolution) can be followed. Variations of this field 
equation are found in many other areas of physics, including electro- 
dynamics and general relativity. The concept of gravitational potential 
energy is presented. The importance of the potential/density relation is 
stressed in the chapter. Simple cases of potential theory are explored: 
point mass, leading to solar system dynamics; the isothermal sphere, 
often used to describe galaxies; the homogeneous sphere, which gives 
a good representation of the central region of our Galaxy; and more 
generic cases frequently used to describe the density profile of galax- 
ies and dark matter halos. The projection of a volume density into 
an (observed) two-dimensional distribution is discussed, along with its 
inverse. 


2.1 Gravitational force 


The dominant force in the evolution of galaxies is gravitation. Over the 
lengthscales and densities typically found in galaxies, a Newtonian ap- 
proach is well justified. Only in chapter 7 will we consider general relativ- 
ity to follow the growth of density fluctuations in the framework of an 
expanding Universe. In its simplest form, the gravitational force on a mass 
mı, exerted by a second mass, mz, separated by a radius vector r12 (with 
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origin at mı) can be written 


mım2- 
37 "12, (2.1) 
M2 


Fi2=-G 


where G is the gravitational constant. This book deals with the role of this 
fundamental force in the formation and evolution of galaxies. We simplify 
the analysis by neglecting other mechanisms such as gas-related phys- 
ics, energy loss due to electromagnetic interations or effects of two-body 
collisions. The latter will be justified in the next chapter. 

Since forces from different sources can be added linearly, we can 
write the net gravitational force on mass m, exerted by a distribution of 


masses {m;} as 
M = Eor, (2.2) 


=-Gmı Ie 
where the second expression takes the system to the limit of a continuous 
distribution, defined by a mass density p(r). An important simplifica- 
tion can be invoked if the distribution of mass is spherically symmetric: 
p(T) = p(r). If the test particle, mı, is located at a distance r from the 
centre of symmetry, we find that all the contributions to the net force from 
the mass outside r cancel out, leading to 


— w [ p(s)4ms2d (2.3) 
r 


Therefore, a test particle moving in a circular orbit with radius r around 
such a distribution of matter will have an orbital speed of 


vn- me 4 p(s)s?ds. (2.4) 


Although stars in galaxies have a wide range of orbital shapes, this simple 
scaling relation is a powerful tool to give us an order of magnitude estim- 
ate, relating speed, mass and size. For instance, the Sun moves with an or- 
bital speed ~220 kms~! around the Galactic centre, located ~8 kpc away 
(~2.5 x 102° m). Therefore, an order of magnitude guess for the Galactic 
mass is Mgix ~ 9.1 x 10!°Mo (where Mo = 2 x 10°° kg is the mass of the 
Sun), which is in the right ballpark. Hereafter, it may help to remem- 
ber the value of the gravitational constant, G, in units more suitable for 
the analysis of galactic systems: G = 4.3 x 1076 kpc(km/s)?/Mpo. It is also 
convenient to know that a speed of 1 km/s is approximately 1 pc/Myr. 
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2.2 The Kepler problem 


The simplest case of a gravitational system is the Kepler problem, con- 
sisting of a star with mass ms at position řs and a planet with mass mp at 
position řp. The forces acting on star and planet are, respectively: 


MsMp 


mi, = GL Fy -7,) 
Irs —Tp| 
(2.5) 
roan cea See 
Mplp = = = 
PR Ifs — fpl? aed: 


Hereafter, we follow the standard notation using a dot (two dots) over 
a quantity to denote its first (second) time derivative (e.g. x= d?x/dt?). 
If we subtract one equation from the other, and define the relative vec- 
tor separation as f = fp — fs, we simplify the two-body problem to a single 
gravitating system with mass M = ms + mp: 

en ee => 1a a =constant = £, (2.6) 

r 2 r 

where the second part is obtained by multiplying both sides by 7 and 
then integrating with respect to time. This first integral represents en- 
ergy conservation. A second conserved quantity can be obtained from the 
definition of angular momentum (per unit mass): 


> 5 


i= xT > 0=Fxr+PFxr=O0, (2.7) 


which cancels for any central force (i.e., as long as the acceleration is 
aligned with the radius). Therefore 2 is a constant vector, and there is no 
loss of generality if we assume the whole trajectory to be confined on the 
XY plane, so that the total angular momentum is: £= £ê,, with £ a con- 
stant. Noting that £ is the area of the parallelogram spanned by ř and ř, 
we arrive at Kepler’s second law: The radius vector joining the star with the 
planet sweeps out equal areas during equal intervals of time. The area swept 
per unit time is dA/dt = ¢/2. Using cylindrical coordinates (R, 6, z) allows 
us to write the energy (per unit mass) as 


R? (2 GM 
E = 2 + Der (R), where PeR) = OR2 = R` (2.8) 


Note this is an equivalent one-dimensional problem involving the radial 
distance from star to planet. 
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Exercise 2.1 


Show that the orbit in a Keplerian potential can be obtained by in- 
tegrating equation 2.8 with respect to the azimuthal angle 0, and 
defining u = 1/R, so that R = —£du /d9, leading to 


fis. ta 
1+ecosé 


which corresponds to a conic section with eccentricity e: 


2E b2 
e=V1+ ape 


and a = £? /GM is the semilatus rectum, directly related to the size 
of the orbit, æ = b? /a, where a and b are the semi-major and semi- 
minor axes of the conic section, respectively. 


Exercise 2.1 confirms Kepler’s first law — The orbit of the planet is a conic 
section with the star at one of the two foci — adding to it how the orbital 
parameters are related to the mass of the star and the energy and angular 
momentum of the orbit. Finally, for a closed orbit, we can derive its period 
either by integrating the energy equation (2.8) with respect to time, or by 
using Kepler’s second law, noting that the area of an ellipse is mab (see 
exercise 2.1), and that therefore the period is: t = 27 ab/¢ = 27a,/aa/¢, 
where substituting a with respect to mass and angular momentum gives 


GM = 02a", (2.9) 


where w= 2z/t is the angular frequency. Equation 2.9 is the mathe- 
matical expression of Kepler’s third law:' The square of the orbital period 
of the planet is proportional to the cube of the semi-major axis of its orbit. 
Moreover, we directly relate this period-orbital size relation to the mass 
of the system. 

We can use the effective potential, defined in equation 2.8, to il- 
lustrate the expected shape of the orbits. The Kepler problem produces 
orbits that are conic sections. The type of orbit depends on the initial 
conditions — that is, the position and velocity given to the planet ‘in the be- 
ginning’. We can use an energy-based argument to determine these orbits. 
Note that equation 2.8 corresponds to a one-dimensional case, namely, 
following the radial distance from the planet to the star. In addition, the 
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Figure 2.1 Kepler’s orbits and the gravitational potential. The bottom 
panel shows the effective potential for a choice of angular momentum 
(£? = 2GMRo/5), and three possible types of orbits defined by the total 
energy: E1, E2 and E3. The orbit for each case is shown in the top three 
panels. 


conserved angular momentum gives us the azimuthal component, via: 
6 =v, /R=/R?. Figure 2.1 shows the effective potential for a specific 
choice of angular momentum, and three possible cases for the total en- 
ergy: E1, E2, E3. Noting that the kinetic energy term is positive definite, we 
find that at the minimum (E1) there is only one possible radial distance, 
so that the planet describes a circular orbit with constant speed (and zero 
eccentricity, i.e., E1 = —G?M*/2¢7). In the region E1 < E < 0, the planet 
is also bound, this time between an aphelion and a perihelion, tracing an 
elliptical orbit. If E > 0, the orbit corresponds to a hyperbola, and at the 
limiting case E = 0, the planet traces a parabola (with eccentricity e = 1). 


Exercise 2.2 


Find the equivalent of Kepler’s three laws of motion if the gravita- 
tional force were to behave like a simple harmonic oscillator; i.e., 
substitute equation 2.1 with 


Fy2 =—Gmjmor12. 
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2.3 Potential theory 


The gravitational force is a vector field: at each point in space, the ac- 
celeration on a test particle is uniquely described by a vector. The grav- 
itational potential is a scalar function that provides the same amount of 
information as the vector field, and it is defined as 


o@=-c ||| pr) ey (2.10) 


IF —F| 


One can show that at each point in space, the gradient of this scalar func- 
tion produces a vector field equivalent to the gravitational force. By using 


the expression 
= 1 -rT 
v| = |==, 2.11 
(r) rai 


we can write the gravitational force on a test particle with mass m as 


ano -Vo P). (2.12) 
m 


Moreover, we can derive a field equation that relates the gravitational 
potential to the mass density. Taking the divergence of equation 2.12: 


Vz- (*) =¢ I “z. (as) p@ ar’, (2.13) 
m i lr’ —r|3 


and using the following identity: 


¥5-( -rT J 3 ition: uma (2.14) 


FP FF 


we note that this expression vanishes when |?” — t| 4 0, so that any contri- 
bution to the integral must come from an infinitesimal volume around r’ = 
T. Therefore, we restrict the volume of integration to a small sphere at” = 
T with radius h > 0. Within the infinitesimal volume of integration, we 
can swap fr and?’ in the gradient operator, allowing us to apply the diver- 
gence theorem, and changing the volume integral into a surface integral: 


> => _ =? . 2È 
Vz- () =—Gp(r) C-D- PS (2.15) 


P-ren FF 
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On the sphere, the surface element is related to the solid angle element 
by d25 = (7 —7)hd2Q, and therefore 


Vz- () =~ 6p) fp P2=-4r60®. (2.16) 


If we write the divergence of the force with respect to the scalar potential, 
we arrive at the gravitational field equation: 


VZO Ë) =4rGo 0), (2.17) 


also called ‘Poissons equation’. Note this equation provides a direct link 
between the source of the field (the density distribution) and the grav- 
itational forces, described by ®(r). This is the Newtonian equivalent 
of Einstein’s general relativistic equation - G,,, =82GT,,, — linking the 
field (through the geometric tensor, Guy) to the mass-energy distibution 
(through the stress-energy tensor, T»). 


Newton's theorems 


Two theorems proved by Newton for spherically symmetric mass distribu- 
tions provide a simple way to solve for the gravitational potential. The first 
theorem states that the gravitational force on a test mass inside a spher- 
ical shell with constant surface density is zero. The second theorem proves 
that the force on a test mass outside this shell is the same as if the whole 
shell mass were reduced to a point mass at the centre. 

There are several ways to solve the theorems. Let us consider a 
test point mass P, inside the shell. Figure 2.2 shows how to solve it via 
solid angles. We compute the gravitational force exerted by the masses 
within two opposite and infinitesimal regions of the shell, which sub- 
tend the same solid angle dQ. When È is defined as the (constant) 
surface mass density of the shell, the mass from each of these two op- 
posed regions is dm, = Ur7dQ/cos6@ and dmz = Erżd2/ cos 0. Note the 
angles are the same, as the triangle formed by the centres of these two 
mass elements and the origin is isosceles. Hence, the forces per unit 
mass are 

F 1 dm, Ed 


m r2 ~ cos 


(2.18) 
Fo _ dmz _ Ed 
m rż ~ cos 
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Figure 2.2 Solution to Newton’s theorem, using a simple geometric 
argument. Note that at an arbitrary point, represented by P, the solid 
angles subtended by dm, and dmz are equal. 


Since they have the same magnitude and opposite direction, they will 
cancel out. Integrating throughout all solid angles, we obtain a net zero 
force. 

We can solve the second theorem by invoking Gauss’s law, mak- 
ing the assumption that, because of the spherical symmetry of the mass 
distribution, the potential (or the magnitude of the force) can depend 
only on the radial distance r. Choosing a spherical surface S, at dis- 
tance r from the origin, enclosing a volume, V, that includes the shell, 
we find 


fp Fy-dé= fff ¥-Fmav=— fff vomav= 


aS(ro) V(ro) V(ro) 


=—4rG I p(r)\dV = —4r GM, 


V (ro) 


(2.19) 


where we also apply Poisson’s equation (2.17). Since the force at S is the 
same in magnitude, we obtain F (r) = GM/r?, as expected if all the mass of 
the shell were located at the origin. 
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Exercise 2.3 


Show that Newton’s first theorem can be solved by integrating 
throughout the spherical shell. An alternative method involves 
showing that the gravitational potential inside any point inside the 
shell is constant. 


2.4 Gravitational potential energy 


We have seen in section 2.3 that the gravitational scalar potential gives a 
complete description of the force vector field, via F=—V®. Such a result 
implies that the work involved in moving a test particle along a closed 
orbit (C) should be zero: 


work= f FaF= || © x -d3=- f| xFo)-dB=0, (2.20) 
C 5 SIS 


where we have applied Stokes’s theorem to the line integral and used the 
fact that the curl of a gradient must vanish. This property characterizes 
conservative forces and allows us to derive the gravitational potential en- 
ergy. Let us assume we have a point mass m at the centre of coordinates, 
and we bring another mass m from infinity, assumed to be initially at rest. 
The net change in total energy will be the work done to bring that particle 
close to the mass at r= 0, namely: 


r 2 2 
work= f F-a¢= | z iac OAR (2.21) 
aT r 


where we assume the particle’s final position at radial distance r from the 
origin. This net change of energy can be associated with an increase in 
the gravitational potential energy: AE = Work. If we do the same with 
an ensemble of masses, we have to add the pairwise contributions to the 
energy change, making sure we do not count twice the same pair, namely: 


w- > Y mof) (2.22) 
= = m T: A 
2 rj 24 S 


where we follow the standard notation using W as the gravitational 
potential energy, and define rj as the separation between the ith and the 
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jth mass. The second part of the equation takes advantage of the definition 
of the gravitational potential seen by the ith particle. This expression can 
be extended to a continuous distribution of mass by replacing the sum 
with a volume integral, using densities instead of masses: 


was I PDPA r. (2.23) 


Based on dimensional grounds, it is often useful to write this expression 
with respect to the total mass of the system (M) and some estimate of its 
size (R): 

GM? 


W=- i> (2.24) 


where the so-called fudge factor (y) is just a numerical value, typically 
not very different from 1. 


Exercise 2.4 


Show that the gravitational potential energy of a sphere with mass 
M, radius R, and constant density is 


Assuming that the solar luminosity Lo ~ 4 x 1033 ergs~! is caused 
by a slow contraction of the star releasing gravitational potential 
energy into heat, derive the expected lifetime of the Sun, and 
contrast your result with its current age ~4.5 Gyr. 


2.5 Potential/density pairs: A few fundamental cases 


Poisson’s equation (equation 2.17) implies that for any mass distribution, 
we can define a gravitational potential from which the force on a test mass 
can be obtained. It is useful to define below a set of density-potential pairs 
that will be very relevant to galaxy dynamics. We assume here spherical 
symmetry. Moreover, it is also practical to derive the radial dependence of 
the orbital speed of a test mass in a circular orbit. One of the key observ- 
ables in galaxies and stellar clusters is the velocity field, from which we 
can derive a large amount of information about the underlying potential. 
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The circular speed gives us a good order of magnitude estimate, and in 
systems where the kinetic energy is dominated by rotation — such as disc 
galaxies — it allows us to probe the mass distribution. 


Point mass: This is the Kepler problem described above, for which the 
density-potential pair is 


GM 
p(r) =MB(r) > &)=-—, (2.25) 


from which we derive an orbital speed v? (r) = GM/r, which corresponds 
to the characteristic orbital speeds of planets in the solar system (where 
M ~ Mo). This is also the typical velocity profile of matter (gas and stars) 
moving close to a supermassive black hole, as found in the centres of 
galaxies. A Keplerian profile is a telltale signature of a highly concentrated 
distribution of mass. 


Homogeneous sphere: The density/potential pair for a sphere with 
radius a and constant density is 


1 
—27 Gpo (a? — at r<a 
(r) = poH(a —r) <> (r) = (2.26) 
P(r) = poH(a —r) (r) Pee 
—_—_—_ r>a. 
3r 


The orbital speed inside the mass distribution rises linearly with radius 
v=r./4xGpo/3, switching over to a Keplerian profile at r > a. Inside the 
sphere, the potential is equivalent to a three-dimensional harmonic oscil- 
lator, with the same spring constant along the three directions, and an 
associated period: t = ./32/Gpo. Note that the dynamical timescale of a 
gravitational system is tayn « /1/Gp. 


Isothermal sphere: An isothermal sphere is a density distribution often 
encountered in gravitating systems. We will see in chapter 3 that it corres- 
ponds to the galaxy dynamics equivalent of a system in thermodynamic 
equilibrium (although there are substantial differences, as we will see in 
that chapter). The density/potential pair is 


P= po (2) <> 9 (r) = 4rGpoa? In(r/a), (2.27) 


and the orbital speed is v2 = 4r Gpoa?, i.e., constant. Both the potential 
and the total mass diverge at large radii. A more realistic version of this 
case involves a truncated profile, where the density vanishes at r > a, and 
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the total mass is therefore M = 47 a", with the potential outside the 
sphere being Keplerian: ®(r > a) = —GM/r. 


Exercise 2.5 


A slight variation of the isothermal sphere is the Jaffe model, where 
the gravitational potential is given by 


x 
B(r) = In | — 
(r) =47Gpoa‘* In ( I =) ; 


and x =r/a. Show that the density of this model is 


p(r) BAET 


and the total mass is M = 47 pọa?. 


NFW profile: The first computer simulations of galaxy formation ex- 
plored the growth of structure under a gravitational potential in an ex- 
panding Universe. As we will see in chapter 7, small density fluctuations 
grow under gravity, creating stable (virialized) structures called ‘halos’. 
The numerical simulations appeared to settle into a universal density 
profile, termed ‘NFW’ after the authors of the seminal work.* An NFW 
distribution is defined by a total (virial) radius (rv), a scale radius (a), 
and a central density po, giving the following density/potential pair: 


po 47 Gpoa? 
= (x)= In(1 2.2 
mE (x) = n(1+ x), (2.28) 


p(x) 


where the dimensionless radial coordinate: x =r/a. The circular speed in 
this potential is 


2 
vx) = an [ha +x)- sl i (2.29) 


The logarithmic slope y =d(ln p)/d lnr of the NFW profile goes from —1 
at r — 0 to —3 at r— œ. At the scale length r =a, the slope reaches the 
value of the isothermal slope, y = —2. Note that since the density scales 
like p xr? at large radii, the cumulative mass diverges logarithmically. 
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Therefore, we also need to truncate the distribution. The outer radius is 
commonly defined as the virial radius, i.e., the region within which all 
particles are bound in virial equilibrium (see section 3.9). The dimen- 
sionless concentration parameter (c) is defined as the ratio between the 
virial radius and the scale radius. This is one of the most studied dens- 
ity distributions, given its relevance in the study of dark matter halos. An 
alternative definition includes free parameters for the behaviour at small 
and large radii in the so-called generalized NFW profile: 


po 


sn F (2.30) 


p(x) 
where (a, )= (1,2) reduces to the NFW profile, (2,0) gives the iso- 
thermal sphere, and (2, 2) is the Jaffe model presented in exercise 2.5. 


Hernquist profile: A similar density distribution was proposed to explain 
the mass distribution in elliptical galaxies.’ In this case, we have 


27r Gpoa? 
a eee i (2.31) 
x(1 +x)’ 1+x 


p(x) 


The two-dimensional projection of this density profile follows quite 
closely the standard de Vaucouleurs (R14) surface brightness profile, as- 
suming a constant M/L ratio. The logarithmic slope changes from —1 to 
—4; therefore the total mass of the distribution is finite (M = 27r pọa°), as 
the density decreases faster than the volume at large radii. The circular 


speed is 
TA cee (2.32) 
v(x) = ,/ — . : 
a 1+x 


Plummer sphere: This potential is especially important in numerical sim- 


ulations of dynamical systems, where the parameter à (softening length) 
eliminates the problem of numerical divergence when computing the po- 
tential for small values of r. In this case, à defines the resolution limit in 
the calculation of gravitational forces. The density/potential pair is 


p(n= (a (1 + r i <> ğ(r)= aM (2.33) 
Ar i3 2 VTZ 42 


Note that as A —> 0, the Plummer sphere resembles a Keplerian profile. At 
small radii (r « A), the model mimics the inside of a sphere with constant 
density. 
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King profile: This model is often adopted to describe the mass distribu- 
tion of stellar clusters, elliptical galaxies and even galaxy clusters. It gives 
a good fit to a numerical method that assumes a Maxwellian velocity dis- 
tribution along with a simple prescription to take into account the ejection 
of stars moving quickly enough to become unbound (see section 8.4): 


„N2 —3/2 
p(r) = po p+(2) l <> 


| 


2.34 
(r/Te) = 


(r) = Do + 4r Gpor? Í = 


Yukawa potential: The Yukawa potential is defined as the Keplerian case 
with an exponential cutoff over a lengthscale ro, namely: 


GM GM 
pr)=-5 e T/T) «> O(r) = en, (2.35) 
ror r 


The Poisson equation for this case gives a negative density! This merely 
states that such a mass distribution that results in exponentially decay- 
ing forces at large distances cannot be accommodated in gravitation; i.e., 
the force of gravity cannot be screened. We will see in chapter 3 that 
this is one of the main reasons why we can consider the motion of a star 
within a large ensemble as produced by the gravitational potential pro- 
duced by the general distribution of stars. The Yukawa potential gives a 
valid description of weak interactions (as in those felt by neutrinos‘). 

Notice that in many of the potentials presented here, the density pro- 
file diverges at small radii. This does not cause a divergence in the amount 
of mass in the central region, as long as the rate of increase of the dens- 
ity is slower than the decrease of volume as r — 0. However, divergent 
central densities (cusps) may be erased by other mechanisms, such as 
the infall of smaller structures (mergers), creating a flatter distribution: a 
core. The core-cusp dilemma in galaxy halos is one of the open problems 
today. 


Exercise 2.6 


Find the circular velocity profile of the last three cases (Plummer, 
King, Yukawa). 
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Figure 2.3 Projection of a spherical distribution onto a two- 
dimensional (observer) plane. 


2.6 Two-dimensional projection 


A comparison of a density distribution with an observed galaxy or cluster 
involves the projection of the three-dimensional volume mass density into 
a two-dimensional surface mass density (or a deprojection if going from 
the observational to the theoretical data). To transform mass into light 
(needed when contrasting models with observational data), we need to 
include the mass to light ratio, Y =M/L. One typically assumes a simple 
scaling of Y with respect to the galactocentric radius. Let us consider here 
the simple case of a spherical galaxy with constant density p (r) = po inside 
r < ro, vanishing outside. We need to integrate along the line of sight (fig- 
ure 2.3) to write the 2D surface brightness, È (R), as a function of the 3D 
density, p(r): 


to ro 2rp(r) 20 
E(R)= yo? t= dr= ,/r2 — R2. 
o I. DA= le e ty (2.36) 


The last step makes the approximation of a constant mass-to-light ratio 
throughout the galaxy. The inverse process, i.e., the deprojection from a 
2D surface brightness profile to a 3D mass density distribution, can also be 
obtained by considering that equation 2.36 describes an Abel transform 
that can be inverted as follows (assuming a constant T): 


Yo [(@ dx dR 
r dR /R2— r2 


p(r)= (2.37) 
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Of course, in practice, this expression may not be practical if the data 
are too noisy or the assumption of some level of symmetry is unreal- 
istic. One useful pair of 2D/3D distributions corresponds to the modified 
Hubble law: 


Saal 9-3/2 
ER) _ 14 (4 o OL 1+ (2 , (2.38) 
Xo ro po ro 


where Yo = 2poro/ Yo. Note that at large radii, the surface brightness de- 
creases as x R~* and the mass density scales like œ r~’. In either case, the 
integrated quantity diverges logarithmically. 


Exercise 2.7 


Consider a galaxy described by a spherically symmetric Hernquist 
profile (equation 2.31), and assume a constant mass to light ra- 
tio, Y. Show that the two-dimensional projected surface brightness 
of this galaxy is 


M (y?+2)F(y)-3 
~ Y2ra? = (1-y?)2” 


2o) 
where y = R/a, and R is the projected radial coordinate, and 


1 
—— sech Yy, forO<y<1 


Ji 
1 


5 sec” ly, forl<y<oo. 
y*-—1 


Notes 


1 Aptly called Kepler’s ‘1-2-3’ law in Misner, Thorne & Wheeler, 1973, Gravitation, 
W. H. Freeman. 

2 Navarro, Frenk & White, 1997 ApJ, 490, 493. 

3 Hernquist, 1990, ApJ, 356, 359. 

4 E.g., Bettini, 2008, Introduction to elementary particle physics, Cambridge. 
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A statistical treatment of 
stellar systems 


A many-body gravitational problem — involving millions or even billions 
of stars — cannot be solved in the same way as the standard two-body 
problem. A statistical treatment is required, following an analogy between 
thermodynamics and statistical physics. The concept of phase space and 
the distribution function is fundamental to a statistical description of a dy- 
namical system. This treatment is valid for stellar systems because of the 
nature of gravity and the characteristic densities in these systems, which 
lead to a long relaxation time, a concept that will be explored in some de- 
tail here. The main equation that governs the evolution of the distribution 
function in a system made up of stars (or dark matter particles) is the colli- 
sionless Boltzmann equation (CBE). Simplifications to the CBE arise from 
symmetries of the system (similarly to Noether’s theorem in theoretical 
physics), leading to the concept of isolating integrals and Jeans theorem 
and allowing us to describe the distribution function in terms of a reduced 
set of variables (the isolating integrals themselves). Simple cases are dis- 
cussed with spherical and axial symmetry. A more pragmatic approach to 
the CBE involves taking moments of velocity, leading to Jeans equations, 
a fundamental workhorse in stellar dynamics. A brief glimpse of the con- 
sequences of perturbatively removing the collisionless behaviour of stars 
is presented, with dynamical friction being one of the most representative 
effects. 


3.1 Phase space 


Since the evolution of a dynamical system is governed by a second order 
differential equation of the position with respect to time (i.e., F=mr), 
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there are two boundary conditions, say at some reference time (to): on 
the position 7(tg), and on the velocity v(to) = 7(to). Therefore, the evol- 
ution can be traced in phase space, comprised of position and velocity 
(or linear momentum). Phase space is therefore six dimensional. Let us 
consider an example of one-dimensional motion subject to a simple har- 
monic potential: ®(x) = akx?. The solution to the equation of motion of 
a particle of mass m in this potential is 


x(t) =A sin(wt + po) | 
; (3.1) 


x(t) = Aw cos(wt + ġo) 


where the angular frequency is w? =k/m; A is the amplitude of the os- 
cillations; and ¢o is a phase term that depends on the initial position and 
velocity of the mass. The trajectory in the corresponding two-dimensional 
phase space is an ellipse. Each ellipse has a unique value of the (con- 
served) energy, E = mA?w?, and the trajectories corresponding to differ- 
ent values of the energy do not cross, creating a nested set of ellipses that 
cover all of phase space. This is an important concept that will be exploited 
when we discuss Jeans theorem below. 


Exercise 3.1 


Find the trajectories in phase space of closed orbits in the Kepler 
problem (see section 2.2). Simplify the 4-dimensional parameter 
space — we can neglect motion along the z-axis — by plotting sub- 
spaces (R,R) and (6,6). 


3.2 The distribution function 


Now we change from a single star moving under a potential, to a many- 
body system comprised of multiple stars. One can study the evolution of 
all the stars in phase space by creating a complex set of trajectories. We 
define the distribution function -f (t, V; t) —as the number of stars located 
within some differential volume element in phase space (dt = d?rd*v), 
at some time, t (hereafter, we use DF to refer to the distribution func- 
tion). After normalization, the distribution function can be considered a 
probability distribution in phase space. For instance, we define the spatial 
number density of stars as 
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v(T; t) =f fE, V; t)dv. (3.2) 
Moreover, we can define observables such as the average velocity as 
ya 1 eee 
(V(r; t) = —— | f(r, v; Dydy, (3.3) 
v(r;t) 
or the velocity dispersion (via the variance) as 


o? T; t) = (WF; 0) — VE; t)? 
ON 1 A SE E E aE . (BR 
= vE; t) | fC, v; t)v d°v VW; t)) (vr; t)). (3.4) 


These will be useful quantities that can be compared with observations. 
For instance, the projection of the average velocity can be extracted from 
spectroscopic observations of galaxies, where the spectral lines appear 
differentially offset depending on their bulk motion, due to the com- 
bination of Doppler-shifted lines for every star in the system. Likewise, 
the velocity dispersion can be inferred from the widths of the spectral 
lines. Most of the time, we assume a steady state, where the distribution 
functions and all the averaged quantities do not have any explicit time 
dependence. 


3.3 Relaxation time 


There is a direct analogy between the motion of stars in a galaxy and the 
motion of particles in a fluid, say, e.g., molecules of nitrogen in a room. 
Although it is possible to envisage the trajectory of a molecule in phase 
space, the number of collisions (understood as interactions with other 
molecules that affect its dynamical state) is so numerous that after a very 
short timescale, these interactions have completely altered its initial dy- 
namical state. This process, called ‘relaxation’, affects all molecules and 
erases all information about their initial dynamical state. This is the main 
reason why it is relatively simple to deal with a complex system made up of 
a large number of particles as a macroscopic (thermodynamic) ensemble. 
In a relaxed system, energy equipartition is achieved, and, in terms of 
the central limit theorem, the distribution function of energies can be de- 
scribed as a Gaussian distribution, with a simple parameter that controls 
its width: 


my2 
fE, V) xe KT. (3.5) 
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Figure 3.1 Computing the relaxation time: simplified (rectilinear) 
trajectory of a star-star encounter. 


This is the Maxwell-Boltzmann distribution. The single parameter 
that uniquely defines this system, temperature, enables us to wholly 
describe its properties. Is it possible to extend this treatment to the dis- 
tribution of stars moving in a galaxy? To answer that question, we need 
to consider the relaxation timescale in typical gravitating systems. 

We can address this problem by studying the trajectory of a star (S1) 
moving in the galaxy. We begin with a single encounter with another star, 
S2 (figure. 3.1). Both stars have mass m. The point of closest approach, 
b, is termed the ‘impact parameter’. We assume here that the deflection 
angle is small, so that the trajectory is approximately a straight line. The 
gravitational force of S2 on Sı is decomposed into longitudinal and tan- 
gential components. There is a net change in the tangential component of 
the velocity from the force: 


-3/2 
Gm2 b 2 t\? 
ee ce he] (3.6) 


b2+x2r b2 b 


This net force will result in a change of the tangential component of 


velocity: 
-3/2 
Gm te vt\? 2Gm 
maar p+(Z) l d= -y (3.7) 


This change in velocity can be interpreted as the force at the point of 
closest distance (r= b) multiplied by the ‘duration’ of the passage (b/v), 
i.e., an impulse. The symmetry of the process results in a cancellation 
of the longitudinal component during the passage (ôv =0). Now we 
assume there are a number of similar passages as Sı moves through- 
out the galaxy. Let us consider the galaxy as a flat system, with N stars 
within a radius R. After a crossing time — defined simply as a typical 
timescale for a star to traverse the length of the galaxy — the number of 
encounters with impact parameter b within an infinitesimal interval db — 
can be written 
N 


én= mpz” bdb. (3.8) 
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The collisions are random events, so that, on average, there is no net 
change in the velocity ((ôv1 ) =0). However, the interactions will impart 
a level of fluctuations, as expected from a random walk process, with a 
root mean square given by 


2Gm\? 2N 
ôv? = (=) pz ab. (3.9) 


The net effect will be obtained by integrating across all impact parameters. 
However, as is always the case in these scattering processes, the contribu- 
tion diverges as b + 0. We simply substitute the lower limit by some Dyin, 
at which our simple approximation of a linear trajectory breaks down. 
bmin can thus be defined as the value of b for which the change in velo- 
city is as large as the actual speed of the incoming star: ôv ı ~ v => bmin = 
Gm/v2. We also change the upper limit to the extent of the galaxy and get 


R Gm\? R 
pe 2p = 
ai= |" wjar—an($") n( =). (3.10) 


Note that the effect of the integration limits appear as a (weakly vary- 
ing) logarithmic trend, usually called the ‘Coulomb logarithm’ (In A). 
Just to get a feeling for the effect of this term, for typical galaxy scales 
(R~ 10 kpc; m=Mo; v ~ 200 km/s) we get bmin ~ 0.02 AU. Interactions 
within this value of b are very unlikely unless the system has density 
~1 (Does which is too high in all stellar systems. Therefore the approx- 
imation is well justified. 

Finally, in order to assess the onset of relaxation, we would need to 
have this level of fluctuation in velocity comparable with a typical velocity 
in the galaxy, roughly v2 = GM/R=GNm/R, where M is the mass of the 
galaxy. Therefore, when relaxation is achieved, after n,e ‘crossings’ of the 
galaxy by the star, 


Av, \2 8ln A 
nel (=) = Nre] SE (3.11) 
v N 


A= =— v? ~N 3.12 
bmin Gm” á £ l 
leading to 
N 
~——., 3.13 
rel 8InN ( ) 
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Table 3.1 Relaxation times of typical stellar systems. 


R v tcross trel 
System N pe km/s yr yr Age/ trel 
Open cluster 100 2 0.5 4x10 107 <1 
Globular cluster 10° 4 10 4x10? 4x108 <10 
Massive galaxy 10/4 104 300 3x10’ 107 1077 
Dwarf galaxy 10? 10 50 2x107 1014 1074 
Galaxy cluster 108 10 108 10° 2x101? 0.1 


This expression can be translated to a timescale by use of the crossing 


time: 
N R 


8InN v’ 
Table 3.1 shows the relaxation times in several stellar systems. In most 
cases, the relaxation time is comparable or even greater than the cur- 
rent age of the Universe. Only in the densest regions (e.g., the cores 
of globular clusters) could one expect the effects from relaxation to 
be significant. In fact, the inherent instability of gravitational forces 
(unscreened, as ‘negative mass’ does not exist!) leads to a rapid col- 
lapse if thermodynamic equilibrium is to be achieved. This issue will 
be explored in more detail in chapter 8. However, for the rest of the 
stellar systems, especially galaxies, we can assume relaxation is never 
achieved. This allows us to adopt the approximation that as a star moves 
in phase space, it is affected only by the (long-range) gravitational poten- 
tial of the system, and not by local interactions. Such behaviour repres- 
ents one of the fundamental properties of dynamical systems moving under 
gravitation. 


(3.14) 


trel = Nrel teross = 


3.4 Local and distant encounters 


We can gain more insight into the role that local and distant encoun- 
ters play in relaxation by giving a more realistic description of a two- 
body interaction. In section 2.2 we derived the trajectories in a Kepler 
potential, finding them to be conic sections. We will focus here on scat- 
tering events, tracing hyperbolic trajectories as illustrated in figure 3.2 
for the interaction between two stars with masses m and M. D is the im- 
pact parameter — the closest approach had there been no gravitational 
attraction — and v» is the original velocity of incoming star m (we will use 
the rest frame of star S2). As r > oo, the azimuthal angle ¢ > 7/24 6/2. 
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Figure 3.2 Computing the relaxation time: encounter of two stars on 
an unbound hyperbolic trajectory. 


Hence 


1+ecos(z/2+4/2)=0 => sin@/2) =~ (3.15) 


z 
Using this expression along with the derivation of the angle 6/2 from the 
triangle XFC, we get 


cot(9/2) = 2 


F (3.16) 
a 


Equations 3.15 and 3.16 can also be used to find the eccentricity: 


_ EN (Y _ [2 (3.17) 
a +(2) ~ V1—cosé’ 


When D —> a, 0 — 7/2. Another useful equation here is the vis viva integ- 
ral found in the Kepler problem. Along the path of the orbit (which, in 
general, could correspond to any conic section), the quantity 


v*=G(M+m) (2+3) (3.18) 
r a 


is conserved (use the negative sign for an elliptical orbit, and the positive 
sign for a hyperbola). In the case r > oo, 


2 _GM+m) Dv? 


Veg a => cot(@/2) — G(M+m) . (3.19) 
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One interesting result from this expression is the concept of ‘gravitational 
focusing’. Let us compare the cross section (i.e., the area within which a 
strong gravitational interaction is likely) with the geometric area of the 
stars (i.e., are and nr, where rı and rg are the radii of stars Sı and S2, 
respectively). A physical collision will take place if the true point of closest 
approach is smaller than the sum of the stellar radii, i.e., a(e — 1) < rı + 
r2. From (3.18) 


eS) i+ 2a = Tmin(24 + Tmin) (3.20) 
= = 7 x 


Fmi i 
min rin 


leads to 
VMAX _ D 


Voo ry +r2 


(3.21) 


For typical values in the solar neighbourhood (vax ~ Vese ~ 620 km s7}; 
væ ~ 30 kms~!) we get a difference between the actual cross section of 
interaction and the simple geometric cross section of 


2 
( Z ) ~ 400, (3.22) 
rı +12 


leading to the idea of focusing: two stars will get closer because of their 
mutual gravitational interaction, increasing their chances of collision by 
a factor ~ 400 with respect to a simple geometric argument (i.e., no 
long-range forces). Taking the local stellar mass density in the solar neigh- 
bourhood as reference, pọ ~ 0.1Mọopc7?, we obtain a typical separation 
between stars à ~ 2 pc. Given that A ~D, we find that rmm ~ 2/20 pe > 
rı +12. Therefore, even if focusing is taken into account, the probability 
of physical collisions is negligible. Only in the densest systems, such as the 
cores of globular clusters, we can expect this type of effect to be relevant. 

Let us now revisit the issue of relaxation time (section 3.3). We could 
adopt the more rigorous treatment of hyperbolic orbits to deal with the 
two-body interaction, and follow a similar treatment as in our rectilin- 
ear case above, integrating with respect to the parameters of the collision 
to create an averaged version of the process. The method is a bit more 
involved, but follows the same argument, estimating a parameter that 
traces the interactions as a stochastic process that introduces a random 
walk in the evolution of the kinetic energy of a test star. We refer the 
reader to the standard treatment of Chandrasekhar,! and show the final 
result for the relaxation time: 


1 /3 (v2)3/2 (3.23) 
TREL = 7 ; 
REL 16 V x nG2m? In (Rv2/GM) 
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where v and m are the typical velocities and masses of stars in the sys- 
tem, and n gives their number density. R is the upper limit to the impact 
parameter, similar to the derivation followed for the rectilinear case, for 
which we assume this term to be roughly comparable to the size of the 
whole system. Just to give a simple estimate, the solar neighbourhood 
(m~1Mo,v~30km sl, n~0.1 pe, R~ 300 pc) has a relaxation time 
TREL ~~ 5 x 1013 yr, therefore much longer than the Hubble time. This 
treatment allows us to determine the contribution to the potential relaxa- 
tion process from nearby or distant ‘collisions’: here we define a ‘collision’ 
as a gravitational interaction that will produce some deflection of the tra- 
jectory of the test star. From equation 3.19 and noting that væ represents 
the typical velocity of a test star, we can write the Coulomb logarithm as 
In(R/a), where a is defined by the hyperbolic trajectory of the star. We can 
now assess the contribution to the relaxation time from collisions with an 
impact parameter D by calculating the ratio 


In(D/a) 


a In(R/a)’ 


(3.24) 
where D/a is related to the deflection angle 0, via equation 3.16. Let 
us consider an example at the highest stellar densities, say, a globular 
cluster core (a~9 x 107° pc; R=Prcore ~ 0.3 pc). In this scenario, colli- 
sions with D/a ~ 1.5 lead to 6 ~ 67° (i.e., very large deflections), and the 
contribution is just X ~ 5%. At weaker deflection angles, say D/a~ 10 
(corresponding to 0 ~ 10°), the contribution increases to X ~ 30%, and 
for very weak scattering events, i.e., those caused by far away stars (say 
D/a~ 500 — 6~ 0.2°), the contribution is X ~ 75%. Therefore, we con- 
clude that in galaxy dynamics, most of the ‘collisions’ between stars that 
eventually cause relaxation are distant. 


3.5 Collisionless Boltzmann equation 


Given that in most circumstances a stellar system can be treated as an 
ensemble of particles without any collisions, with individual trajectories 
corresponding to the underlying gravitational potential, we can describe 
this motion in phase space according to a simple differential equation. 
Consider the motion of these stars in phase space as some sort of ‘fluid’, 
moving through a differential element of phase space, dt = d?rd?v. We 
will use Cartesian coordinates for simplicity, so that the differential ele- 
ment is a six-dimensional hypercube. We want to describe the change in 
phase space density, f (f, V; t), following two steps. First we write down 
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the flux of stars through the opposite faces in position of the hypercube. 
For instance, through the faces normal to the êy direction, we have a flux 
difference of 


an 


which, combining all three spatial directions, leads to the first part of the 
variation of the distribution function: 


ne au || ay dedtudv aw = — ar, (3.25) 


— V;(fv)dt = —¥- Vef dt. (3.26) 


The second part of the variation of the distribution function arises from 
the flux through the ‘velocity’ faces of the hypercube. For instance, the 
contribution along the é, direction (i.e., the x-component of the velo- 


city) is 
0 
bi sie Cap tu |} aydedudvay 8 CD ge = 


dt dt ə 


_ af a® 


ae re —dt. (3.27) 


The final expression uses the fact that the x-component of the acceleration 
can be written with respect to the gradient of the gravitational potential 
(from chapter 2). Note that this step also needs to assume the potential 
does not depend on the velocity (always true in Newtonian gravity). The 
variation of the number of stars within this volume of phase space should 
equal the implicit dependence of the distribution function with time, 9; f. 
The combination of all these terms results in a continuity equation: 


Fs -Vf — VO@)- Vif =0. (3.28) 


The above equation is the collisionless Boltzmann equation (hereafter 
CBE). We can understand this equation if we consider that in fluid dynam- 
ics, the Lagrangian derivative is defined as the rate of change of a physical 
quantity in the fluid as one follows its trajectory. In an ordinary fluid, this 
derivative, applied to the fluid density (pọ) can be written 


Do =%p+¥-Vp=0. (3.29) 


We can extend this concept to phase space, so that the CBE can be written 
with respect to a Lagrangian derivative as 


Df=0. (3.30) 
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Therefore, the motion of stars in phase space is such that as we follow 
the trajectory of one of these stars, we do not find any local variation of 
the stellar number density (Liouville’s theorem). An alternative way to 
arrive at this interpretation is by looking at the CBE as an equation in par- 
tial derivatives in seven coordinates (phase space plus time), leading to 
six subsidiary, first order, ordinary differential equations: 


dt dx dy dz dua Ww 
oe ee) ee) ee 


These correspond to the standard equations of motion of a test particle 
under the gravitational potential of the system. 


(3.31) 


Exercise 3.2 


Give the two main reasons why the distribution of stars in the Milky 
Way galaxy cannot be described by analogy with respect to the dis- 
tribution of air molecules in a room. How about the description of 
stars in a dense globular cluster? 


3.6 Isolating integrals: Jeans theorem 


The CBE governs the behaviour of the distribution function in phase 
space. It encodes all information about the dynamical system, although 
normally we use the velocity moments to make comparisons between 
models and observations. It is not trivial to determine the distribution 
function that corresponds to a stellar system, and the use of conserved 
quantities is fundamental to be able to operate with this formalism. In 
this context, we can separate conserved quantities in two categories: isol- 
ating and nonisolating integrals of motion. The best way to visualize the 
difference between them is to take pairs of dynamical variables (i.e., posi- 
tion and its associated momentum) as action-angle pairs. It is easiest to 
picture the evolution of such a pair: a trajectory will cover the surface 
of a torus (see figure 3.3). The trajectory of a star with a given con- 
stant of motion may be as in panel (a) of the figure, where the orbit is 
closed and the whole trajectory does not cover the whole surface of the 
torus. Another option is shown in panel (b), where the trajectory does 
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Figure 3.3 Illustration of the difference between nonisolating (a) and 
isolating (b) integrals. 


not close and ergodically covers the whole surface. Keeping in mind that 
no two orbits should cross (as they correspond to different integrals of 
motion), one could say that by choosing a conserved quantity that be- 
haves as in panel (b), we can separate this part of phase space volume 
from the rest, by choosing that specific value of the integral. Different 
values of the integral will define different tori. A choice of the related 
constant of motion can thus be used to ‘compartmentalize’ phase space. 
These constants of motion are called ‘isolating integrals’. In the more gen- 
eral case, a given isolating integral will divide phase space into disjoint 
hypersurfaces, a second isolating integral will define a different set of hy- 
persurfaces, and so on, so that one could fully describe the distribution 
function with respect to these integrals, rather than the standard way 
using position and velocity (Jeans theorem). Moreover, these isolating 
integrals emanate from a symmetry of the system, allowing us to produce 
educated guesses for the distribution function of a stellar system, as we 
will see in section 3.7. 

The isolating integrals most often found in stellar systems are en- 
ergy, expected in stationary systems, i.e., those without an explicit time 
dependence in the distribution function; total angular momentum (and 
its components) if the system is spherically symmetric; or a single com- 
ponent of angular momentum if only rotational symmetry is present 
about one of the axes. For instance, in the Milky Way galaxy, to low- 
est approximation, we may consider the system to be in steady state and 
with cylindrical symmetry, so that E (energy) and J, (the projection of 
angular momentum along the axis of rotation) are isolating integrals. In 
chapter 4 we will see that the motion of stars in the Galaxy suggests an ad- 
ditional third integral. Although finding isolating integrals is non trivial, 
one can make use of the properties of the system to produce educated 
guesses. 
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3.7 Examples of distribution functions 


We will show here the properties of some distribution functions (DFs) fre- 
quently used to describe stellar systems. All systems are assumed to be in 
a stationary state. Hence, energy (E) is always an isolating integral. Note 
that if E is the only isolating integral, the velocity distribution must be iso- 
tropic. This stems from the fact that E=v*/2+ ® (P), so that the second 
order moment of the j-th component of velocity is 


ve) = Ta FE v-F(E)d°y, (3.32) 


leading to (v2(7)) = (vy (r)) = v2 ®©). Therefore, if the observed velocity 
distribution is not isotropic, there must be an additional isolating integ- 
ral that breaks this symmetry. Such is the case with motion in an axially 
symmetric system (e.g., a rotating disc), where the projection of angular 
momentum along the symmetry axis (say Jz =rv,) breaks the symmetry 
with respect to v2. 

When defining distribution functions, one commonly encounters an 
energy floor, corresponding to the minimum of the effective potential 
(®o), from which a relative potential can be defined as 


UF) =—O(7) + Go, (3.33) 


as well as a relative energy: 


- v2 - 
EM=—Z +O). (3.34) 


Polytrope 


The polytrope model defines the distribution function as some power law 
of the relative energy: 


ee! 
Reale" > a (3.35) 
0, E <0, 


where n is the polytropic index. The density can be written with respect 
to the relative potential by integrating in velocity space: 


Sau y2 
p=m ||| rEv=4am | F(Y -7V v= =Y", (3.36) 
0 
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where the coefficient cn is defined as 


(n— 3)! 


c= mA(2z)?/2 7 


(3.37) 
Note that for this coefficient to converge, we need n>1/2. There- 
fore the case n=O (i.e., a homogeneous polytrope) is not physically 
possible. We can write the equation that leads to the radial depend- 
ence of the density by using Poisson’s equations in spherical polar 
coordinates: 

dY 


V?Y = —4r Gp => == (« =) + 4 Gc," =0, (3.38) 
r2 dr dr 


which is the Lane-Emden equation, typically found in the analysis of stel- 
lar structure, when describing a system in hydrostatic equilibrium, as 
illustrated in exercise 3.3. 


Exercise 3.3 


By assuming hydrostatic equilibrium in a spherical distribution of 
gas with a polytropic equation of state p x p”, show that we recover 
the equivalent case to a collisionless polytrope (i.e., equation 3.38), 
with y=1+1/n. 


Two interesting cases arise for the polytrope model: n=5 produces the 
Plummer sphere presented in chapter 2, and n-> oo leads to the iso- 
thermal model, presented below. 


Isothermal sphere 


In stellar dynamics, the isothermal sphere is the thermodynamic equi- 
valent to a Maxwell-Boltzmann distribution. Replacing the kinetic en- 
ergy by the total energy, and the temperature (which behaves as an 
energy scale, kT) by a velocity dispersion parameter (øo), the isothermal 
sphere is 


E 
Aes? , E>0 
F(é)= (3.39) 


E<0. 


© 
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Exercise 3.4 


Show that for the isothermal sphere distribution function defined 
in equation 3.39, the density profile is p (r) « 1/r?. 


The isothermal sphere is one of the most useful distributions in 
gravitational dynamics. Note that the velocity of a circular orbit is con- 
stant (vo): 


GM(<r) G [fT K v2 

2 2 2 0 

v = atl 4rs as 4nGK =v% = p(r) een 
(3.40) 


Also note that the cumulative mass profile increases linearly with radius. 
Therefore, in order to avoid a divergent mass, the density is truncated at 
some radius ro. This defines the so-called Truncated Singular Isothermal 
Sphere. In spiral galaxies, the velocity profile outside of the central 
regions is roughly flat, suggesting a quasi-isothermal distribution. 


Mestel disc 


Mestel disc? is a distribution related to the dynamics in disc galaxies 
such as our Milky Way. This distribution features a net bulk rotation, and 
the stars move within a very thin disc. Therefore, the three-dimensional 
mass density, o(r), is replaced by a two-dimensional density, £ (R). Note 
that, hereafter, we distinguish between the 3D radial distance (r) and 
the 2D projection (R). The equivalent of an isothermal sphere in a 2D 
system is given by È (R) « 1/R. Denoting vo as the radius-independent ve- 
locity of any circular orbit in this system, we can write the surface mass 
density as 


zR = (3.41) 
anh 
and the associated relative potential as 
2, [R 
Y(R)=—vln( — }, (3.42) 
Ro 


where we assume the gravitational potential, ®, vanishes at R= Ro. The 
cumulative mass profile also increases linearly with radius and calls for a 
truncated model, where the mass density vanishes at R > Ro. 

Following Jeans theorem, we assume that energy (E) and the 
vertical component of angular momentum (Jz) are isolating integrals 
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(we assume the disc is located on the xy horizontal plane). We follow the 
previous ansatz of an isothermal distribution for the energy and add a 
power law term with respect to angular momentum: 


Te NE fe. 
a( =) eo, Jes DAL=0 
0, Jg<OVE<O. 


(3.43) 


The condition on the sign of Jz eliminates counterrotating stars. From the 
DF, we can derive the surface mass density: 


E(R)= I F(E,J;)d*v = 
love) +00 Rv, J 

=A dv dv), | —— X 

f = [. l (= p 


oO 
a | 
as 
ol 
a| + 
N| < 
HN 
Qa j< 
blob 
5 
fo 
D| > 
Se 
Ld 


leading to 


1. / Ro \1(/R\ 0" 
ER) =Ao2n2227(4 3.45 
(R) =Ao“n Cr) Fave Ro (3.45) 


This result is compatible with our original estimate of the surface mass 
density (equation 3.41) if 


q= 2-1. (3.46) 


This expression allows us to consider two different regimes: if q is a large 
number, the system is a cold dynamical disc, where the kinetic energy is 
mostly in bulk rotational motion (vo > o). If q > —1, then we have a hot 
dynamical disc, and the kinetic energy budget is dominated by random 
motion (loosely called ‘pressure’). 


Osipkov-Merritt models 


An extension that includes angular momentum in the DF to a spherical 
system is given by the Osipkov-Merritt models.* We consider a 3D system 
where the isolating integrals are energy (£) and angular momentum (J), 
but the dependence of the DF is given by a variable Q: 

J2 


fE, Vv) =F(Q), where Q=€+ =>. (3.47) 
2r 


a 
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Note that this expression is motivated by the contribution of angular 
momentum to the energy. For instance, in the Kepler problem (see 
chapter 2), the total energy can be written 

R GM P 


E=— 3 : 
7 R + SRZ (3.48) 


In this case, ra is a free parameter that relates to the anisotropy of the 
orbits, as we will see below. We separate the velocity vector into a radial 
and a tangential component: v = v,ê, + Vt. The mass density is 


p= || FQdy=20 f av, f vrr @= 


2 0 pe * 
a. | dQ/2Q— FQ, 
1+% 4* 


(3.49) 


where ®(r) is the gravitational potential. We are interested in the an- 
isotropy of the velocity distribution, so let us compute the radial and 
tangential components of the velocity dispersion tensor: 


1 2w f° 
pvp) == Sy f dQ [2(Q— $)]”? FQ), 
3145 Jo 
ra 
(3.50) 
2 2m f° 
a ae f dQ [2(Q— ©)]*/? F(Q). 
1+ z ® 
Therefore, the level of anisotropy is given by the ‘beta’ parameter: 
2 2 
pate eS (3.51) 


2 y2 27 
Oo; rera 


In this model, the central regions (r « ra) have an isotropic velocity dis- 
tribution (£ ~ 0), whereas the outskirts (r >> ra) feature radial anisotropy 
(8 ~ 1). This behaviour is typical of galaxies, where the outer regions 
have a more significant presence of stars on orbits which are very close 
to radial trajectories, reflecting events in the mass assembly history. Re- 
member that in galaxy dynamics, the relaxation times are very long; 
hence stellar orbits retain their past dynamical history for long periods 
of time. 
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3.8 Jeans equations 


The CBE is defined in the six-dimensional phase space, making it hard 
to exploit with observational data, which offer instead an ‘integrated’ 
mapping of this equation. The Jeans equations provide a more tractable 
set of equations for comparisons with data. They consist of moments of 
velocity: 


Jeans equations => i dvv? (CBE), (3.52) 


where n is the moment order. To avoid overloading the expressions in the 
rest of this chapter, we will use a single integral sign to denote ‘velocity- 
volume’ integrals. Let us begin with the zeroth order moment: 


I afdv+ Y j VidxfdPv — Y dx, ® f æ fdv=0. (3.53) 


The last term vanishes: apply the divergence theorem, which reduces the 
volume integral to a surface integral evaluated at v > +00, where the dis- 
tribution function should vanish. Note we have also pulled the spatial and 
time dependence out of the velocity integrals as these are independent 
variables. The remaining integrals refer to the number density v (r; t) and 
average velocity (v;(r;t)), leading to a continuity equation, this time on 
averaged variables: 


dv + > ðe (v(vi)) =0. (3.54) 


We can relate the stellar number density, v, to the mass density in stars, 
p =mv, where m is the mass of a single star. 

The first order moment of the j-th component of velocity is obtained 
in an analogous way: 


f dfvjd?v + >D f vivjðx; fdv — 5 Ox, D [rs =0. (3.55) 
i i 


Applying again the divergence theorem to the last term and assuming a 
vanishing surface integral (in velocity space) leads to the following: 


Ie(V(Yj)) + X dx; (v(ViV;)) + V9; = 0. (3.56) 
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The second term is usually written with respect to the covariance of the 
velocity distribution (the velocity dispersion tensor): 


oF = ((vi — (Vi) Yj — (Vj))) = (Vij) — (Vi) v), (3.57) 
leading to 


VdE(Y;) + VV) » Vx; (Yj) + vay D+ X anvo) =O. (3.58) 


l 


Jeans equations in spherical coordinates 


We will show below the derivation of a very useful representation of Jeans 
equations. We assume here spherical symmetry, and consider only the 
equation created by multiplying the CBE by the radial component of the 
velocity. The three components of velocity and acceleration in spherical 
coordinates are shown in the equation in exercise 3.5. 


Exercise 3.5 


Show that in spherical polar coordinates (r,0,¢), where 6 and ¢ are 
the polar and azimuthal angles, respectively, the three compon- 
ents of velocity are v-=7; Vg = r and Vo = rọ sin, and the three 
components of the acceleration are 

a, =? —r6* — rd? sin? 0 

dg =rÖ + 276 —r¢? sin 0 cos 8 

ag = rọ sin @ + 27¢ sin 0 + 2rd cos 6 


(Hint: find first the time derivatives of the unit vectors: Ê., ĉo, êp). 


We begin with the most general expression for the CBE: 


af+ > rdf +Y viðnf =0, (3.59) 


and use the vector equation Ÿ = —V ® in spherical polar coordinates, not- 
ing that, by the assumed symmetry, the potential depends only on the 
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radial component, leading to 


dd  %+V% 
Vr = — — 

r dr r 
oo v cot 0 — vr vg ; (3.60) 
Vg = —____ 

r 

. Vor + VoVo cot O 
Vo = 


r 


We replace these expressions into the CBE, and follow the same pro- 
cedure as above, multiplying by the radial component of the velocity 
and integrating in all of velocity space. Let us write the expression 
in full: 


a: / wef +3, / BPw?f + dp | dêvv,vef +g / Bw, vef 


2 2 2 
Vo tv d vi cot — vrvo 
+ fom (8-2) aps omg 


tø 
= f dew, eee rere a =O. (3.61) 


Some of the integrals can be readily evaluated as averages of velocities 
(using equations 3.2, 3.3, 3.4). The last three terms require integration 
by parts, where we are left with integrals such as 


I dygivi anf = fgv iS- J Prfyyg(vis---), (3.62) 


and the first term vanishes, as it is evaluated in the limit towards in- 
finite velocities, where we expect the distribution function to decrease 
faster than the increase in any other velocity-related expression — to 
avoid, for instance, an infinite velocity average or velocity dispersion. The 
expression simplifies to 


Ie(v (Vr) + Ae (v (Ve) + 3o © (vrvo)) + Ap V (Vrvg)) 


2) 4 (y2 2 (3.63) 
(2 (v2) Fjot 


DGLP) 
r dr r 


Since the system is spherically symmetric, we expect the average velo- 
cities to vanish ((v;)=0), allowing us to write o? = (v?). The velocity 
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dispersion tensor can have only two non-zero components: Orr =o; 
and ogg = op =0_. The off-diagonal terms vanish. Therefore, the final 
equation is 


1a(vo2) 202 dọ 
red aa ay (3.64) 
v ar r dr 


where we have defined the anisotropy parameter £ =1 — (o? /oŻ). We 
will also see in chapter 4 the equivalent version of Jeans equations in 
cylindrical symmetry. 


Exercise 3.6 


Let us describe a dark matter halo as a spherically symmetric distri- 
bution with a constant and isotropic velocity dispersion o. The halo 
includes a baryon, diffuse, gaseous component without a tempera- 
ture gradient. The dark matter density dominates the mass budget 
(ppm > pg). With the use of Jeans equation, assuming a stationary 
state and no bulk motion, show that the gas density profile and the 
dark matter density profile are related via: 

umpo? 

kT 


pgr) <[ppm(r)]", where 7 = 


If the anisotropy parameter, £, is a function of the radial coordinate, the 
solution to the homogeneous version of equation 3.64 is 


(vo?) =K exp (-2 a E ap) 3 (3.65) 
o P 


and X is a constant that trivially corresponds to the value of (vo,2) at the 
origin of coordinates. The general solution can be obtained by promoting 
the constant to a radial function K(r), leading to 


© | GM(<s)v(s) ° BP) 


where we have applied the boundary condition K(r) — 0 as r— oo. 
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3.9 The virial theorem 


The Jeans equations were obtained by removing information about in- 
dividual velocities and replacing them by averaged quantities. Further- 
more, one can integrate the Jeans equations over the spatial coordinates, 
effectively removing all space-related information and leading to an equa- 
tion that deals with the energy balance in the dynamical system. This is 
the virial theorem: 


virial theorem => l dx, (Jeans equations). (3.67) 


Let us start from the equation corresponding to the j-th component of 
velocity (equation 3.56): 


| sae(otv)yx= = 5 fads (otvivj)) d?x — I PXk Ix, dx= 
' (3.68) 
= > | p(Viv;) din 2x + Wag = 2Kig + Wijs 


where we have integrated by parts the first term on the RHS and ap- 
plied the divergence theorem, so that the surface integral vanishes as 
x— œ. The two integrals on the RHS are symmetric with respect to 
the indices k,j, and we define them as K and W, representing the 
kinetic and potential energy tensors, respectively. The kinetic energy 
tensor can be split into two parts related to bulk motion (7) and bulk 
motion (TT): 


1 


1 
Kj = 5 / plvyj)dex = 5 f p [vo (v) + og] dx= Tij +M. (3.69) 


The gravitational potential energy tensor can be duly identified if we write 


it as 
ə ne i 
Wy =G / p(X)xXk— | L eu dx | dx= 


Ox; Ix —x’| 
(3.70) 
=6 ff osc A aude. 


The expression is symmetric with respect to the integrating variables 
x and x’. We can symmetrize it, i.e., add a new term switching these 
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variables, and multiply all by a factor of +, to get 


= Xk) 0G — 34) 
Wy =—> 5 fog oO rss zB 7 dxd. (3.71) 
The trace of this expression is 
Í 1 
trw= > Wii = 5 f p(x) (x)d2x, (3.72) 


j=l 


which is the standard gravitational potential energy (see equation 2.23). 


Exercise 3.7 


Elliptical galaxies are believed to form from the mergers of pro- 
genitor systems that can be either discs or other ellipticals. Let 
us consider two elliptical galaxies with the same mass (M), size 
(ro) and velocity dispersion (oo), slowly approaching from infinity 
and eventually merging to form a more massive galaxy that settles 
into virial equilibrium. Show that the total energy of the system is 
Eror = —3Mo¢. Find the velocity dispersion, size and density of the 
new galaxy, assuming that its mass distribution is identical to that 
of the original galaxies. 


3.10 Beyond the collisionless Boltzmann equation: 
The Fokker-Planck equation 


When collisions need to be taken into account (e.g., when describing the 
dense, central region of a globular cluster), the CBE is not an adequate 
description of the phase space density, and an additional term needs to 
be considered: 

Df əf 


of 
Ta ore Vf- VD- Vf = (4) (3.73) 


The term on the RHS takes into account the change in phase space dens- 
ity from collisions. We can describe this term by a probability distribution 
function, such that p(v, AY)d AY gives the probability that a particle with 
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velocity v experiences a shift in velocity Av within a differential element 
dAv over a time interval At. This probability gives information about 
individual encounters. According to this definition, we can write 


+00 
fO,t+ A= f@— AY, t)p(v — AV, Av)d? Av, (3.74) 


—00 


thereby transforming the Boltzmann equation into a rather complicated 
integro-differential equation (called the ‘master equation’). 

Even in dense stellar systems, we are still in a regime where the re- 
laxation time is long enough. Therefore Av «v during a crossing time, 
and we can perform a Taylor expansion (LHS in At and RHS in Ay) and 
truncate at the second order; this is the Fokker-Planck approximation: 


af 3 a(F(AVi)) 1 À FF (AvjAY;)) 
2 NEE, 3.75 
( Ja ee ( ) 


ot OVjOV; 


with the first and second order moments of the change in velocity given by 


(Av) = / pv, AV) Avid? Av, 
l ; (3.76) 


(AvyjAyj) = / pv, Av) AyjAvjd? Av, 


commonly known as encounter integrals, or diffusion coefficients. The 
first term in the RHS of equation 3.75 describes a drag process (dynamical 
friction). The second term describes a diffusion process in velocity space 
(compare it with the diffusion equation 3;f = —kVf). Note that the more 
general treatment of the Fokker-Planck equation (3.73 and 3.75) would 
involve the diffusion coefficients to be integrated over velocity and pos- 
ition. However, we apply the local approximation, where the individual 
interactions take place over regions much smaller than the size of the 
system. Hence, we can simplify the treatment of interactions as a super- 
position of Keplerian hyperbole unaffected by the global potential of the 
system. 


Dynamical friction 


In equation 3.75 we came across the dynamical friction term. Physically, 
one can think of it as a drag force (i.e., velocity-dependent) where the 
passage of a mass, say an object like a black hole, a globular cluster or 
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an incoming galaxy, affects the surrounding region, altering the potential 
that, in turn, produces a back reaction. It can be represented pictorially 
as a wake behind the moving mass that changes the local density of stars, 
introducing an additional force. There is a useful expression to describe 
dynamical friction, due to Chandrasekhar. We assume that a mass M 
moves through a homogeneous sea of stars, each of mass m <M. The 
motion can be simplified as a combination of hyperbolic trajectories in 
a Kepler potential (see section 2.2). Section 3.4 describes the details of 
a hyperbolic trajectory, but we neglected there the effect of the scatter- 
ing event on mass M. By considering conservation of momentum, we can 
relate the reaction on M from the scattering of m as imparting a velocity 
change: 


Avy = — (ax) AVm, (3.77) 


and we know from section 3.4 that the effect on m is a rotation of the 
incoming velocity vector Və% by an angle 6. Now we use the frame of ref- 
erence of star M, and add many such scattering events as it traverses the 
field of stars with mass m. For a single scattering event, the net change in 
velocity is 


5 (3.78) 


where we split the contribution into a parallel and a perpendicular com- 
ponent, with respect to the motion of M. The symmetry of the problem 
will give, on average, zero net change in Av], and so we need only to 
consider the parallel contribution. The deflection angle can be written 


2 
cos =1- 5, (3.79) 
e 


where e is the eccentricity of the hyperbolic trajectory. We also introduce 
here the impact parameter, i.e., the distance of closest approach during 
the scattering event -b = a (e — 1) - and use the expression of eccentricity 
from equation 3.17, leading to 


2mv D\? 7 2mv b2y4 - 
faa CO 2 = CO oe) . ; 
Av = = a+(2) l Pe re P z] (3.80) 


Now we can add up all the contributions from each star as an integral 


with respect to the impact parameter b, including the number density of 
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stars in this homogeneous sea, by use of the distribution function, that, 
by definition, can depend only on velocity, f Vm): 


dm _ 


bax 
a Om — 7m) | 2b Av, (b)dbf Om)d? vm. (3.81) 
JO 


In this expression, we have replaced Yæ by the relative velocity between 
mass M and one of the stars, i.e., Vm — Vy, so that we can add up the 
individual contributions. The range of integration in b extends to some 
maximum value byax, normally represented by the size of the stellar dis- 
tribution. Since the integral produces a logarithmic term (the Coulomb 
logarithm), our choice is not critical in the derivation: 
dm 2 Vm — VM - 33 
— = 4r G^ (M + m)m = ln Af (Vm) d? vm, (3.82) 
dt Vm — Vm]? 
where A = byaxVoo/G(M + m), and we have simplified the Coulomb log- 
arithm by assuming In(1 + A?) ~ 21n A. This expression corresponds to a 
specific choice of velocity of the sea of stars. Therefore, we need to integ- 
rate it in Vm. By use of the vector expression encountered in equation 2.11, 
we can write 
3 > Vin = YM ind > 
Ja Vmf Vm) sz = VyyH m), (3.83) 
i vm — Yml? 
and H is called the first Rosenbluth potential. If we assume this potential 
to be isotropic, the gradient can be written: 


VH = (ðv HYm/vm, (3.84) 
arriving at the final expression for the dynamical friction effect, namely: 


dv v rM 
TM —167r?G? (M +m)mln ee f fm) v2 dym. (3.85) 
dt Vm JO 


The integral represents the number density of stars with speeds below vy. 
The force is directed against the motion of mass M and depends on its 
velocity. If we assume M >> m, and a Maxwellian distribution of velocities 
for the sea of stars, with dispersion o, we find 
dým 2 ÝM 
—— = —4r G Mp InA-,-2(5); (3.86) 
dt vM 
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where p is the mass density of stars, s=vy/o V2, and the function g(s) 
can be written with respect to the error function: 


2s _ 2 
EEn ae ; (3.87) 


Dynamical friction is responsible, for instance, for the orbital decay of 
globular clusters as they move around the galaxy. Note that this expres- 
sion assumes a point mass M moving through a homogeneous sea of stars 
with mass m, but a similar result can be derived for more realistic scen- 
arios, such as the merging of two galaxies. From equation 3.86 we expect 
a drag force proportional to M? and inversely proportional to vå. There- 
fore, dynamical friction is more efficient when the incoming galaxy in a 
merger is more massive: major mergers — those where the masses of both 
galaxies are comparable — proceed faster than minor mergers, but fast 
encounters — those in high-density environments — will imply a weaker 
contribution from dynamical friction. 


Notes 


1 Chandrasekhar, 1960, Principles of stellar dynamics, Dover, p. 48. 

2 Throughout this book, unless explicitely noted, these isolating integrals are defined 
per unit mass. 

3 This is a simplified version of the model presented in Mestel, 1963, MNRAS, 126, 
553. 

4 Merritt, 1985, AJ, 90, 1027. 

5 Chandrasekhar, 1943, ApJ, 97, 255. 
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4 
Understanding our Galaxy 


The Milky Way Galaxy can be described, to the lowest order, as a system 
of stars in a thin, differentially rotating disc. A more detailed description 
includes additional components: bulge, thick disc, stellar halo and dark 
matter halo, with the latter dominating the overall gravitational poten- 
tial. Each component is in a different dynamical state, imprinted during 
the formation history of the Galaxy. After a brief introduction to the main 
properties of our Galaxy, including basic details of positional astronomy 
(galactic coordinates and the local standard of rest), the chapter dives into 
the standard approach of measuring galactic rotation and to the definition 
of Oort’s constants. A perturbation analysis of idealized circular orbits 
leads to harmonic oscillations, and to the concept of epicyclic motion. In 
addition to rotation, stars experience vertical motion, perpendicular to 
the plane of the disc. This motion is also treated in a perturbative way. In 
the context of the distribution function, one can identify the separation 
of motion on the plane and perpendicular to it as a fundamental isolating 
integral (the so-called third integral). The chapter finishes with applica- 
tions of Jeans equations to the Milky Way and introduces typical functions 
adopted to describe the potential of the Galaxy. 


4.1 General description of the Galaxy 


Our Galaxy (see figure 4.1) - commonly written with a capital ‘G’, to dif- 
ferentiate it from other galaxies — is a late-type system, i.e., a disc galaxy, 
with a fairly typical luminosity, around 2 x 10!°Lo (where Lo represents 
the solar luminosity). The chemistry and kinematics of its stellar pop- 
ulations reflect several components, with a different formation history, 
namely: 


Bulge/Bar: Stars in the central region of the Galaxy have different 
properties than those in the solar neighbourhood, in terms of the 
distribution of their ages and chemical composition, as well as their 
kinematics. In contrast to disc stars, bulge stars have a higher frac- 
tion of their kinetic energy in the form of random motion. In addi- 
tion, there is a bar component, with a significant fraction of radial 
motion. We know from chapter 3 that stellar orbits keep track of the 
past dynamical history. Stars in the bulge encode its past formation 
both through the orbits and the chemical composition, revealing an 
early and intense formation process (hence the old ages and high 
metallicity of bulge stars). The stellar mass in the bulge comprises 
roughly one-fourth of the total. 

Thin disc: This is the component where the solar system is em- 
bedded, and it dominates the stellar mass budget of the Galaxy. 
It can be described by a flat structure, with vertical scalelength 
~0.3kpc, and an exponentially decaying profile along the radial 
direction, with scalelength ~ 2.5 kpc, stretching out to 15-20 kpc. 
The Sun is located at a distance R ~ 8 kpc from the galactic centre. 
The motion of stars in the thin disc is dominated by bulk ro- 
tation, with a roughly constant tangential velocity in the neigh- 
bourhood of the Sun, around 220kms~!, therefore orbiting the 
Galaxy with a period of ~220 Myr. Stellar populations of different 
ages feature different velocity dispersions and vertical scalelengths, 
with the youngest stars having smaller vertical motion and lower 
dispersion. 

Thick disc: This is the second disc-like component of the Galaxy, 
with a significantly larger vertical extent (~ 1 kpc), featuring older 
stars with different chemical composition (higher [Mg/Fe], see sec- 
tion 6.4), although the transition between the thin and the thick disc 
is very gradual, so it is difficult to identify stars — especially close 
to the Galactic plane — as belonging to either the thin or the thick 
disc. 

Stellar halo: Made up of preferentially old and metal-poor stars, 
moving with high velocities and populating a large spheroidal 
volume. Its mass contribution to the Galaxy is relatively minor, 
amounting to about 1 per cent of the total mass in stars. 

Dark matter halo: Although the observations can target only 
the stellar or gaseous component, we should keep in mind that 
the whole system is embedded in a much larger structure made 
up of dark matter particles. In chapter 7 we will explore in 
more detail the role of dark matter in galaxy formation. We can 
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Figure 4.1 Gaia/DR2 view of the Milky Way galaxy. 
(Source: courtesy European Space Agency/Gaia/Data Processing and 
Analysis Consortium.) 


simplify this component as a spheroidal or triaxial system that 
extends much farther out that any of the above components, to 
~150 kpc, comprising about 95 per cent of the total mass content 
of the Galaxy. However, note that as it occupies a larger volume, 
the dark matter mass density is rather low. Such low densities 
explain why, for instance, within the comparatively minute volume 
spanned by the solar system, the contribution from dark matter is 
negligible. 


The coordinate system 


The position of stars in the sky is determined by assuming a so-called 
celestial sphere; i.e., irrespective of their distance to us, we assign two 
angular coordinates, equivalent to the longitude-latitude system that pin- 
points a location on Earth. These coordinates are defined in various ways, 
depending on the adopted frame of reference. The most important ones 
in astrophysics are the following: 


1. Alt-Azimuthal coordinates: The simplest choice uses the local 
frame of reference; i.e., the vertical angle above the horizon is the 
elevation, and the horizontal angle - measured from due South — is 
the azimuth. These coordinates are useful to locate a star in the sky 
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at night, but depend on the time of observation and the location of 
the observer. 

2. Equatorial coordinates: The Earth’s rotation determines the frame 
of reference. The horizon is replaced by the celestial equator, 
i.e., the projection of the Earth’s equator on the celestial sphere. 
The ‘vertical’ angle with respect to the equator is the declination. 
The ‘horizontal’ angle, measured along the equator from the ver- 
nal point — the place where the Sun crosses the celestial equator 
at the (Northern) Spring equinox! — is the right ascension. This is 
the standard coordinate system for all extragalactic astrophysical 
sources, as it is independent of the time or location of observation, 
although it is affected by the long-term precession and nutation 
of the Earth’s rotation axis, requiring the definition of a reference 
epoch (e.g., J2000.0). 

3. Ecliptic coordinates: This system uses the ecliptic as the frame of 
reference, i.e., the projection on the celestial sphere of the Earth’s 
orbit around the Sun. The ecliptic latitude and longitude are the 
equivalent angles, where the zero point on the ecliptic is also the 
vernal point (note the ecliptic plane and the celestial equator inter- 
sect at the vernal point, and its opposite position, i.e., the location 
of the Sun six months later, at the (Northern) Autumn equinox). 

4. Galactic coordinates: This is the most convenient way to describe 
the location of sources in our Galaxy, also given as longitude and 
latitude. The reference is the projection of the galactic plane on the 
sky, and the zero point of longitude is the position of the galactic 
centre. 


One can translate the coordinates of a celestial object among these differ- 
ent systems via a set of rotations. In addition to the standard location on 
the sky, the positions of galactic objects are also typically described by a 
right-handed Cartesian {X,Y,Z} coordinate system, where the origin is at 
the solar position, the X axis points towards the Galactic centre, the Y axis 
points along the direction of rotation and the Z axis is perpendicular to 
the disc plane. 


Exercise 4.1 


Where on Earth would the Alt-Azimuthal and equatorial coordinate 
systems be one and the same? 
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4.2 Differential rotation in the Galaxy 


The motion of stars in the solar neighbourhood and, by extension, in 
the disc can be described by a differential rotational motion. Most of the 
kinetic energy in the stars of late-type (i.e., disc) galaxies is in the form 
of bulk rotation. Note this treatment cannot apply to bulge or halo stars. 
Historically, this description of motion provided the first steps towards 
our understanding of the structure of the Galaxy, with astronomers like 
Oort, Kapteyn and Lindblad pioneering the field. Two different types of 
stellar motion can be observed. On the celestial sphere we can measure 
the proper motion with respect to distant, i.e., fixed, sources. This angu- 
lar motion on the celestial sphere can be mapped onto a physical motion 
on the normal plane with respect to our line of sight, hence defining a 
tangential velocity. Along the line of sight it is also possible to measure 
a radial velocity via the Doppler shift of spectral lines. Both components 
need to be referred to a fixed frame of reference. The first correction needs 
to take into account the motion of the Earth around the Sun, which in- 
troduces a seasonal variation along the ecliptic of Avorg = +30 km s71. 
Furthermore, these heliocentric velocities need to take into account the 
orbital motion of the Sun around the Galactic centre. We therefore need 
to define a Galactic local standard of rest. 


Local standard of rest 


In order to understand the motion of stars in the Galaxy, we need to 
describe their observed velocities as measured from a frame considered 
“at rest” with respect to the bulk rotational motion of the Galaxy. There 
are two main definitions of such a local standard of rest (LSR), and it is 
important to note the difference: 


e Kinematic LSR (kLSR): Defined by a frame of reference located at 
the position of the Sun, moving with the average velocity of the stars 
in the neighbourhood. This is the easiest system from the observa- 
tional point of view, as it requires only knowledge of the distribution 
of stellar velocities of nearby stars. 

e Dynamical LSR (dLSR): Defined as a frame of reference also located 
at the solar position, but moving in a circular orbit on the Galactic 
plane, whose velocity is dictated by the gravitational potential. This 
definition is the natural one from a theoretical point of view, but we 
will see that kLSR and dLSR are not the same and that they move 
apart. This effect is termed ‘asymmetric drift’, and we will describe 
it below. 
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Figure 4.2 Outline of the system involving a nearby star (S), the 
Sun (L), representing the observer, and the Galactic centre (G). 


Lowest order motion: Oort’s constants 


The simplest model assumes that all disc stars move on the Galactic 
plane along circular orbits following the gravitational potential. An 
arbitrary potential is adopted, although for simplicity we will con- 
sider spherical symmetry. Note at this point that the potential relates 
to the total mass content. Therefore, the dark matter halo must be 
included. 

We want to describe the motion of a disc star from our vantage 
point. The observed motions will be given with respect to the Sun, as- 
sumed to follow a simple circular orbit around the Galaxy. Figure 4.2 
outlines the measurement, where L is our (heliocentric) reference, S is 
the target star, and G is the Galactic centre. We denote 69 and 6 the 
angular velocity of the circular orbits of Sun and star, respectively. We 
begin by writing the radial and tangential motion of the stars, relative 
to L: 

v, = RÊ cosa — Roĝo sin £, 
: . (4.1) 
vt = RO sina — RoG cos £ 


Eliminating « from segments GP and LP and making a Taylor expansion, 
assuming S is close to L (i.e., d« Ro), we find, to the lowest order: 


y= =- sin 2£, P 
vem GRO Tnt 02 [o+ l)a 
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These two equations are usually written 
v, =Ad sin 2£, 
>, (4.3) 
v: =Ad cos 22 + Bd 


where A and B are Oort’s constants of Galactic rotation, defined at 
R=Ro: 


(4.4) 


1fv dv 
B=-3 | R . ( aR ) 
Ro 
Oort’s ‘constants’ therefore represent a first order description of the rota- 
tion curve of the Galaxy. Note that A and B in general are not constants, 
but depend on the location of the observer — much in the same way as 


Hubble’s constant depends on the cosmological epoch of the observer (see 
chapter 7). 


Exercise 4.2 


Show that an isothermal sphere yields A = —B =1/,Qo; a spherical 
distribution with homogeneous density results in A= 0, B= —Qo; 
and a Keplerian potential gives A =3/,Qo0, B= —1/,Qo, where Qo = 
vo/Ro is the angular speed measured at R= Ro. 


In the solar neighbourhood, Oort’s constants are:? A = 15.3 + 0.4km s71 
kpce~! and B = —11.9 + 0.4km s~! kpc™!, which, according to the three 
options in exercise 4.2, suggests a density distribution closer to an iso- 
thermal sphere. 


Exercise 4.3 


Leaving aside potential issues of stability, contrast the values of 
Oort’s constants for a homogeneous spherical distribution and a 
homogeneous disc (i.e., flat) distribution. 
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A more general approach to Oort’s constants 


The derivation shown above is restricted to the case where the orbits are 
perfectly circular. Hence, at each point the velocity is purely tangential 
and constant at fixed radius. We can extend this simple model to a more 
general case — still assuming motion on the plane — if we write the velocity, 
measured from the dLSR, as 


VR) =H(Ro) - (R— Ro) + O(IR— Rol”). (4.5) 


Note this is a general Taylor expansion of the velocity field, retaining only 
the linear terms. The matrix H can be interpreted as a Jacobian between 
velocity and position. In Cartesian coordinates, 


2 ƏxVx OyVx. K+C A-B 

H(Ro) = = 6 

so a l ve K-c)’ ce) 
0 


where the second expression defines four Oort’s constants, including our 
previous definitions of A and B, this time extended to the general case, 
where a radial component of the velocity is present, and where v} can 
depend on the azimuthal angle. We used the simplified notation x = 0/d,, 
and so on. In cylindrical coordinates, the generalized Oort’s constants are 


aeiy 13 
= Vv Yy 
RYL — gI 


2| R Ro 
1 VL 1 

B= a ð 
3 | R RYL + R ml 


(4.7) 


1 VI 1 
K=- — + ORV) + = 3V L 
Ro 


2 [R R 


Note that A and B are equivalent to our previous definition in equation 4.4, 
when vj =0 and 3v1 = 0. Furthermore, in the simplified case C =K =0. 
If we keep the axial symmetry but introduce a non-zero radial term, e.g., 
radial migration, Oort’s constants A and B remain unchanged, but C and 
K will be non-zero, although potentially small with respect to A and B, 
as long as this migration term is small compared to the rotation velocity. 
In the solar neighbourhood, the quoted values of the new constants? are: 
C=—3.2+0.4kms—! kpc—!; and K = —3.3 + 0.6 km s7} kpc™t. 
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Exercise 4.4 


Consider a galaxy with a homogeneous density distribution, where 
stars are assumed to move on circular orbits. Via the Doppler effect, 
we observe the radial velocity of a star at distance d and Galactic 
longitude (¢). The distance from the observer (O) to the Galactic 
centre (GC) is Ro. Find the radial velocity as a function of dis- 
tance (d) from the observer. Do the same for an isothermal sphere 
distribution. 


The next order: Epicycles 


Let us go back to our original assumption of uniform circular motion (i.e., 
C= K =O), but now introduce small perturbations about this motion. As 
is common practice in physics, a first order correction —i.e., retaining only 
linear terms in the equations of motion — will lead to harmonic oscillations 
that will allow us to understand the more complex trajectories in the real 
system. We still constrain the orbits on the Galactic plane, but the radial 
and tangential components are displaced with respect to the ‘equilibrium’ 
case of a circular orbit with radius Ro (see figure 4.3): 
E=R-— Ro 
| ; (4.8) 


n = Ro (0 — 80) 
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View looking from NGP 


s / 


Figure 4.3 Coordinate system used to derive epicyclic motion. 


Both perturbed coordinates are assumed to be small with respect to Ro, 
allowing us to perform a Taylor expansion and retain only the linear 
terms. The position (Ro, 00, z = 0), also called the ‘guiding centre’, follows 
a circular orbit with angular speed 


dðo Vio 
ae (4.9) 
The acceleration of the motion can be approximated by the centripetal 


acceleration of the circular orbit: 


v 
a~ —êr—=. .10 
a~ -êr (4.10) 


Now, the acceleration in cylindrical coordinates (see exercise 4.5) is 
a= (È — RÖ” )êr + (2RO + RÖ)èg + Bez. (4.11) 


Assuming the dominant component is only the centripetal term leads to 


a ; v4 
R-R6?=-=. (4.12) 
R 
Furthermore, conservation of J, (the vertical projection of the angular 
momentum vector) gives 


R*6 =Rov_ 0. (4.13) 


Equation 4.13 can be used to eliminate the 6 term in equation 4.12. 

Finally, we write the equation in terms of the small displacements (equa- 
tion 4.8) giving 
n. R&? v2 

=2 40 _ 4, (4.14) 
R3 R 


After some algebra, expanding vı (R) in a Taylor series and keeping only 


the first order terms, we can write this equation with respect to Oort’s 
constants of rotation: 


E£ =4B(A —B)£, (4.15) 
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which is the equation of a simple harmonic oscillator, with solution: 
&(t) =H sink (t — to). (4.16) 


This motion includes a new frequency, termed the ‘epicyclic frequency’, 
namely: 
= 2,/—B(A—B). (4.17) 


The rate of change of £ is therefore 


(YPRo 
=. 


é=Hk cos x(t —to) ~H 


(4.18) 


To describe the tangential motion (along the 7 direction), we begin with 
the conservation of angular momentum (equation 4.13), written in terms 


of £, to lowest order: 
R 2 
ġ ~ 0 (1- x): (4.19) 


The first term in brackets simply tracks the dLSR, whereas the second term 
gives the motion relative to the dLSR. Taking only the latter, and writing 
down the solution to &(t) gives: 


2V 1,0 
D2 
Rok 


Ad = (Yi) Ro sin «(t — to). (4.20) 


The tangential separation of the star and the dLSR is then 


. 2V10 
n= f Poadae= Ry (vi Ro cos«(t— to), (4.21) 


which can be written in terms of Oort’s B constant 


(idr 


2B 


n= cos x(t — tg). (4.22) 


Using equation 4.22 along with equations 4.16 and 4.18, we find the re- 
lative motion of nearby stars around the dLSR as a retrograde epicyclic 
orbit, with axial ratio B = An/Aé = —x /2B (see figure 4.4). The axis ratio 


can also be written 
/ A 
= 1 — —. 2 
B B (4.23) 


The dLSR acts as the guiding centre of the epicycle. Notice that the char- 
acteristics of the epicyclic orbits depend purely on the behaviour of the 
Galactic rotation curve at Ro. 
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Figure 4.4 Epicyclic motion. 


Exercise 4.5 

Show that in cylindrical coordinates, the velocity can be written 
y= Rép + ROE, + 2, 

and the acceleration is 


a= (R — RŽ )êr + 2(R — È + RO)eq + Bey. 


Using the derived values of A and B, we find the epicycles have a period 
Tr = 27 /k ~ 170 Myr and axis ratio 6 = 1.51 (ellipse elongated along the 
tangential direction). Given that the ‘Galactic year’ is 226 Myr, nearby 
stars take 0.75 galactic rotation periods to complete an epicycle around 
the dLSR. 

Noting that the typical speed along the radial direction is (vi )Ro 
~ 30km/s, we find the dimensions of the ellipse are 1.26 kpc along the 
n direction, and 0.83 kpc along the & direction. So, most stars in the 
solar neighbourhood come from within a radial distance approximately 
Ro +1 kpc (from the Galactic centre). 

Following a star along its epicycle, we note that the ratio of average 
speeds is 


(ivil) (d&/dt) (vı)o(—2B)(cos«t) 1 
— — = = = 0.68. 4.24 
(vı —vyL,ol)  (dn/dt) @ı)ox (sin xt) B on 
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However, this measurement cannot be derived from the observations. In- 
stead, we compute the average over many stars moving in the solar neigh- 
bourhood, each having its own epicentric radius. 

Some stars will have guiding centres Rg > Ro, i.e., farther away from 
the position of the dLSR. Those stars will move in the solar vicinity faster 
along the tangential direction, because of angular momentum conserva- 
tion. Other stars will have guiding centres Rg < Ro, showing up in the 
solar neighbourhood with slower tangential speeds. The relative tangen- 
tial speed can be written (using conservation of angular momentum, and 
keeping only lowest order terms in £): 


(Rg) 
Rg 


ù = Ro (Ô — ĝo) ~ Ro jito + 2é |== 2B. (4.25) 


Taking an average of the square of the tangential velocity gives 


+2) _ ap2/¢2 £2 22 
(9°) = 4B°(E*) = Pe (°) = a gE ). (4.26) 
Hence, the ratio of RMS speeds is 
(vt) -B 
iy AEB ~ 0.46, (4.27) 


which is compatible with the observations. Note the epicycles are more 
elongated along the tangential direction, whereas the velocity distribu- 
tion is elongated along the radial direction. 

If we note that the number density of stars is higher at small radii, 
we expect more stars with Rg < Ro, i.e., moving with slower tangential 
speeds in the solar neighbourhood: this produces the so-called asym- 
metric drift. This effect is stronger in older stellar populations (they 
have higher random speeds, so their orbits deviate further from circular 
motion). 


Exercise 4.6 


A galaxy is described by a homogeneous, spherical distribution of 
stars with total mass Mo and radius Ro. Find the ratio between the 
angular velocity and the epicyclic frequency. Can such a potential 
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produce closed orbits (i.e., orbits that repeat themselves after a 
number of periods)? 


4.3 Vertical motion 


We can make a simple estimate of the force normal to the Galactic plane by 
assuming an infinitely thin, plane-parallel slab with constant density, p. 
Let us assume the vertical motion is independent of the epicyclic motion 
on the plane (figure 4.5). At point P above the plane, the acceleration due 
to the volume element dAdz is 


- GodzdA, 
dis" ep, (4.28) 
r 
whose s component is 
4 GodAd 
da, = da - ê; =—= 9S = —Gpdzd9, (4.29) 
r r 


and we use the solid angle of the volume element as viewed from P: 
dQ=-—-. (4.30) 


Given that the slab is infinite and homogeneous, we can trivially integrate 
in solid angle, getting the total vertical acceleration from the whole mass 


plane as 
daz = — 2x Gpdz. (4.31) 


Note that the acceleration is independent of z. The disc of the Galaxy has 
a vertical extent, so we consider as the next approximation, two thin slabs 
symmetrically located above and below the z = 0 plane, at some distance 
+z, with the same density: o(—z) = p (2). 


So, for z(P) >z daz = —4r Gpdz, 
(4.32) 
whereas for z(P) < z daz =0, 


and all material farther out from the Galactic plane does not affect the 
acceleration at P. Hence, the total acceleration can be written as 


z=P 
äz= -4c | p(z)dz. (4.33) 
0 
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Figure 4.5 Vertical motion. 


Near z=0 a constant density profile can be adopted (po ~ pg). This is a 
good approximation in the region |z| < 100 pc. Hence, at small separa- 
tions from the Galactic plane, the acceleration will increase linearly with 
slope 47 Gpo. If we define È as the surface mass density of the disc, i.e., 


+2DISC 

L= I p(z)dz, (4.34) 
—ZDISC 

then the acceleration at |z|, far from the disc, will be -27GX. 

However, the observed vertical acceleration is found to increase be- 
cause of an extended component (a halo). By examining the density of 
stars in the halo, it is evident that most of the matter causing this beha- 
viour must be in the form of dark matter. Hence, the density of dark matter 
can also be derived from vertical motions of stars in the Galaxy. The slope 
da,/dz at z= 2 kpc is especially useful. However, if the observed z is very 
large, the infinite plane approximation will break down. 

Let us now derive the vertical motion for small displacements. If z is 
small, one can assume p ~ pọ, and 

dz 


ae = —4r Gpoz, (4.35) 


az= 
describing a simple harmonic oscillator with solution: 


z = zmax sin[v(t — to)]. (4.36) 
Differentiating, equation 4.35 gives 


v? = 4r Gpo frequency, 
on lz (4.37) 
t= — = |- period. 
v Gpo 
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For po ~ 0.15Mopc~?, the period of the vertical motion is tz = 68 Myr. 
We can now estimate the maximum excursion from the plane, zmax, for 
stars in harmonic motion. We just need to evaluate the velocity of the stars 
when crossing the Galactic plane. At z ~ 0, we observe vz ~ 6 km/s. Com- 
paring z(t) and dz/dt we find oe = vZmax, giving zmax ~ 67 pc, reached 
after tg ~ 17 Myr. 

Therefore, stars in our neighbourhood are expected to probe a tor- 
oidal region centred on the LSR, with radial extent AR = +250 pc over a 
timescale t, = 175 Myr; and a vertical extent Az = +67 pc, with timescale 
Tz = 68 Myr. If any of the periods involved in this motion coincide, e.g., 
tr(R)=1/, of a Galactic year, then resonance ensues, leading to density 
anomalies such as spiral density waves (see chapter 5). 


4.4 The collisionless Boltzmann equation in 
galactic coordinates 


In order to apply the CBE to the Milky Way, we take advantage of the 


axisymmetric shape, writing the equation in cylindrical coordinates 
(R, ġ, z), where the distribution function reads 


fE, V =f(R, p, Z, VI, VL, V25 0), (4.38) 
and the CBE can be written 


Of +Rapf + Pagf + zaf + Ha f +01 dy f +V2dy,f =0. (4.39) 


The general equations of motion for a test particle in cylindrical coordin- 
ates — note ¢ (azimuthal angle)¢ ® (potential) — are 


ab v2 
MR TR 
eee MPL g, (4.40) 
Riad R 
b a® 
a ag 
leading to 
YL vi 
Bef + vy nf + dof +Vedef + | E — RD ) af 
(4.41) 


1 
E (vivi +3 P) dy f — 0,04,,f =0. 
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We make the following assumptions: (1) the system is in steady state 
(af = 0); (2) there is full axial symmetry: 0g = 0 (for any function). Now 
the CBE is: 


VV 
R 


2 
v 
vp def + vsdef — (Ane = )anf dy f —dz@dy,f=0. (4.42) 
This CBE (an equation in partial derivatives involving five variables) can 
be broken into four first-order differential (subsidiary) equations: 


dR dz _ dv _ dv, E dvz (4.43) 


i Yz (4 — ae) (-"h) (- 22) 


Therefore, a maximum of four isolating integrals can be invoked. Two are 
obvious from the symmetries of time (steady state) and rotation about the 
z axis, namely, energy and the projection along é, of angular momentum: 


1 
T=E=5 (vi + +v) + OR,2), (4.44) 


Tə = Jz = Rv. (4.45) 


If these two are the only isolating integrals, Jeans theorem implies that 
fR, Z, Vi VL, Vz) = F(E, Jz). (4.46) 


Notice vj and v; enter only into Z; and only as a term vi +7v2); hence they 
should be interchangeable in the analysis, leading to the same velocity 
dispersion: 

(vq) = (V2). (4.47) 
The observations of the motion of stars in the Milky Way lead to a different 


outcome: 
0.4v?) = (v2), (4.48) 


suggesting an additional isolating integral in the system. 


Vertical motion and the third isolating integral 


The description of stellar trajectories in the galaxy suggests it is possible 
to separate between motion on the disc and motion perpendicular to it, 
in the flattened potential of a disc. This behaviour reveals the presence of 
a third isolating integral. Close to the Galactic plane, the z component of 
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the force — which depends on the local mass density — is almost decoupled 
from the R component — which depends on the enclosed mass. We have 
seen this behaviour in sections 4.2 and 4.3, where the epicycles and the 
vertical motion of stars moving close to the Galactic plane can be treated 
as independent harmonic oscillators. The potential is, therefore, expected 
to be separable into an R-dependent function and az-dependent function, 
namely: 

PR, z) ~ P1 (R) + P22), (4.49) 


which allows us to solve one of the subsidiary equations of the CBE (see 
equation 4.43): 


dz dvz 
vz  —d®2/dz’ oe 
resulting in an integral of motion: 
l2 
T3 = z” + ®2(2), (4.51) 


which can be defined as the energy involving vertical motion (being in- 
dependently conserved with respect to the energy related to rotation). If 
we take this as a third isolating integral, the distribution function for a 
flattened axisymmetric distribution (such as our Galaxy) now reads 


fR, Z, Vj VL, Vz) = F(E, Jz, 13), (4.52) 


whose dependence on v, is no longer the same as on vj, breaking the 
symmetry between (vi) and (v2). 


4.5 Application of Jeans equations 


There are several applications of Jeans equations to the Milky Way. In 
chapter 3 we presented the derivation of Jeans equations as a mar- 
ginalization of the CBE in velocity space. The continuity equation (see 
equation 3.54) can be written in cylindrical coordinates as follows: 


av 4 1 a(Rv(yy)) 2 O(V(Vz)) 


ot R OR 0z es 


The derivation of the set of Jeans equations involving the first order mo- 
ment of the velocities (from the general expression in equation 3.58) is 
more involved but rather straightforward, and leads to the following: 
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am), 8OO) | BOOM) (2 =, a) e 


at IR dz R OR 
aw) wvv) IVY) 2v 
w ta a "Rre 
Iw) aww Iw (vivz) 3P 
I aR a +( R ~)=0 
(4.54) 


We will show below two applications of Jeans equations to the dynamics 
of the Milky Way. 


Vertical motion and Jeans equation 


The first example concerns the equation corresponding to the vertical 
component of velocity. Let us assume stationary state — i.e., neglecting 
the partial derivative in time — and decorrelation between motion on the 
plane and vertical motion — leading to (vv) = 0. The latter approxima- 
tion is especially relevant for motion close to the Galactic plane and is well 
justified by our definition of the third integral. The last equation in 4.54 
now reads 


1 a(v(02)) ad dz 


v az əz d? 


anc | p(s)ds, (4.55) 
0 


where we have used our results from section 4.3 about the vertical acceler- 
ation in the infinite plane-parallel slab approximation. In this expression, 
the density and RMS of the vertical velocity can be observed, allowing us 
to determine the local mass density. Oort obtained a value of the dens- 
ity p(Ro,2=0) ~ 0.15 Mo pe~* (Oort limit). Notice the density requires 
a double differentiation, whereas the surface mass density needs only 
one derivative (i.e., it is less uncertain). Oort’s estimate was (Rog, |z| < 
0.7 kpc) ~ 90 Mo pe~?. 


Asymmetric drift 


The kLSR and dLSR (see section 4.2) are expected to be one and the 
same if all stars move along circular orbits dictated by the gravitational 
potential. However, their positions slowly drift with time if we adopt epi- 
cyclic motion instead. We can give an adequate representation of this 
effect. Let us consider the radial velocity component of Jeans equation 
(the first equation in 4.54), with the additional assumption of symmetry 
of the density profile along the vertical direction. For positions close to 
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the plane, we have zv ~ 0, and so 


RƏM alive) 2 2, _3ð 


We define the azimuthal velocity dispersion as 
of = (v1 - v1)?). (4.57) 


After some algebra, and noting that the velocity of the dLSR is v2 = 
R(d®)/dR, we arrive at the asymmetric drift equation: 


3 Inwo?)) ( “| R “cou 
J í l 


(4.58) 
Note that v, tracks the motion of the dLSR, whereas (v; } corresponds to 
the average circular velocity for a number of stars (i.e., the KLSR). Let us 
define vig sp = (V1) and vaLsr = Ve. Also, since |vasr — VkLsR| & VaLSR, We 
can write 


võLsR — Vise ~ 2VALSR (VALSR — VKLSR), (4.59) 
and so j 
(v7) 
VdLSR — VkLSR ~ T J eek IF (4.60) 


The terms within the brackets are the same ones as in equation 4.58. The 
first term in brackets is the dominant one. The last term can be neglected 
for our purposes. As the density decreases outwards, the first term in 
brackets is positive, i.e., V4LsR > VkLsr, and so, the kLSR lags behind the 
dLSR. Empirically, vatsp — VkLSR ~ (vi) /D, with D=120kms7!. 


4.6 The potential of the Galaxy 


The observed rotation velocity profile of disc galaxies (see figure 4.6) 
can be described to the lowest order by three regimes, as sketched in 
the right panel of the figure. In the central region, we have a linearly 
rising portion, which is the standard rigid body motion expected when 
the mass density is constant (see section 2.5). Most of the rotation curve 
of the galaxy is relatively flat, suggesting an isothermal density profile. 
The sketch also shows the expected Keplerian decrease of the rotation 
velocity in the regions where the mass density drops to zero. However, 
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Figure 4.6 Comparison of the observed rotation velocity profile of 
a disc galaxy (left, NGC3198), with a simple model involving three 
regions (right). The data for NGC3198 also show the individual 
contribution of gas, stars and dark matter to the rotation curve. 
(Source: adapted from Dutton et al., 2005, ApJ, 619, 218.) 


observations of dynamical tracers in the outer regions of galaxies (us- 
ing, e.g., halo stars, planetary nebula or globular clusters) do not seem 
to reach such regions, suggesting that the dark matter halo is much more 
extended than any of these tracers. The observational data in the figure 
show the individual contribution of gas and stars to the rotation curve. 
Note that if we write the circular velocity as vi = GM(< R)/R, we can add 
the individual contributions to the orbital velocity in quadrature, namely: 
vI =Vetars + Vgas + Ym: Aside from such simple models to describe the 
galaxy, one can create more complex mass distributions to mimic the 
rotation curve (and any other dynamical information, such as Oort’s con- 
stants) as accurately as possible. The Schmidt model represents one of the 
early attempts,* with a linear superposition of nonhomogeneous spher- 
oidal mass distributions, along with a point mass, leading to a rotation 
profile: 


v2 (R) = EE + an om a fat 2z (4.61) 


where each component is defined by a spheroid with eccentricity e; and 
density p;(a), and a is the semi-major axis of each spheroid. In fact, this 
model built upon a previous one defined by Oort that assumed only con- 
stant density spheroids.” More complex models are currently adopted to 
describe our Galaxy. As an example, we show below the potential-density 
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Figure 4.7 Miyamoto-Nagai model of the Galaxy. The left panel shows 
the best-fit rotation curve that consists of two components, with para- 
meters a, b and M shown in equation 4.62: one corresponding to the 
bulge (subindex 1) and another one representing the galactic disc 
(subindex 2). The panel on the right shows a vertical cut of the density 
distribution, with contours in units of Mọ pe~*. (Source: Miyamoto & 
Nagai, 1985, PASJ, 27, 533.) 


pair corresponding to the Miyamoto-Nagai model:° 


GM 
O(R,z) = : 
Je + (a2 + vB? +2) 
(4.62) 
ses b?M aR? + [a+ 3V 22 + b?][a + Vz? +b?" 


4m {R24 [a+ J 22 + b2]2}5/2(z2 + b2)3/2 


Note the limit a > 0 represents a spherical Plummer model, and b — 0 de- 
scribes avery flattened potential. For a specific combination of two models 
(i.e. two choices of a, b and the associated mass), one can obtain a good 
fit to the rotation curve of our Galaxy (see figure 4.7). 


Notes 


1 I.e., when the Sun moves from the Southern to the Northern celestial hemisphere. 
2 Bovy, 2017, MNRAS, 468, 63. 

3 Bovy, 2017, MNRAS, 468, 63. 

4 Schmidt, 1956, Bulletin of the Astronomical Institutes of the Netherlands, 13, 15. 
5 Oort, 1952, ApJ, 116, 233. 

6 Miyamoto & Nagai, 1985, PASJ, 27, 533. 
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5 


Specific aspects of disc and 
elliptical galaxies 


This chapter presents a few properties of the dynamics of stars in the 
two general morphological classes of galaxies: disc and elliptical systems. 
Their different dynamical states (discs supported by rotation; ellipticals 
supported mostly by random motion) and surface brightness radial gradi- 
ents (with discs having exponential distributions and ellipticals featuring 
steeper profiles) are telltale signatures of their different formation histor- 
ies. We begin with a general overview of their main properties. We explore 
the scaling relations in galaxies (Tully-Fisher in discs, Fundamental Plane 
and its projections in ellipticals), and their connection with the virial the- 
orem, as well as other physical processes. Regarding ellipticals, we will 
briefly present the concept of random motion support and its connection 
with the observed ellipticities of galaxies. A brief section is included on 
spiral structure in disc galaxies, revisiting the epicyclic motion of stars 
on the disc, and presenting a scenario (density wave theory), where or- 
bital resonance can lead to the rigid body rotation of a pattern (the spiral 
structure). 


5.1 ‘Hot’ versus ‘Cold’ dynamical systems 


We saw in chapter 1 that galaxies can be morphologically classified into 
three main families: spheroidal, disc-like and irregular. The morpholo- 
gical appearance of a galaxy gives a direct link to its dynamical state. 
Detailed studies of the kinematics of the gaseous and stellar compon- 
ents reveal a substantial difference. If we split the total kinetic energy of 
the stars in a galaxy between ordered motion (rotation) and random mo- 
tion (velocity dispersion), we find that spheroidal galaxies are supported 


mostly by velocity dispersion, whereas disc galaxies keep most of their 
kinetic energy in the form of bulk rotation. If we use a simple analogy 
between velocity dispersion and temperature, one can say that this clas- 
sification divides the sample into ‘hot’ and ‘cold’ dynamical systems, or 
slow and fast rotators, respectively. There are two parameters often used 
to quantify this split: the ratio (v/o), typically measured within an effect- 
ive radius (i.e., the region within which half of the total flux is observed), 
and the specific angular momentum parameter: 


(R\v|) 
(RVVv2 +02) 


A= (5.1) 


Incidentally, the traditional morphological separation between elliptic- 
als and spirals as presented by the Hubble tuning fork diagram does not 
convey the true dynamical evolution of galaxies, and an alternative classi- 
fication scheme was proposed to take into account angular momentum in 
addition to morphology.! Note, for instance, that only one-third of nearby 
ellipticals shows a hot dynamical state. 


Late-type galaxies 


Late-type galaxies (LTGs) consist of the family of disc galaxies, featuring 
spiral arms of any type (termed Sa, Sb, Sc, Sd in the classical notation) and 
barred systems (SBa, SBb, SBc, SBd). The subclass label from a to d refers 
to a gradation from close, tightly wound spiral arms (a) to more open 
arms (d). Lenticular galaxies (SO), although classified as disc systems, 
are not considered late-type galaxies. In addition to the disc component, 
late-types feature a central bulge. The so-called classic bulges have dy- 
namical and population properties similar to those of ellipticals. However, 
a second type of bulge (pseudo-bulge) is expected to form instead, after 
the onset of dynamical instabilities within the galaxy. The bulge-to-disc 
ratio — which can be measured either by mass or luminosity — decreases 
along the Hubble sequence (largest for SO, smallest for Sd galaxies). 

If we invoke conservation of angular momentum, the collapse and 
subsequent cooling of the baryonic material is expected to produce a 
galaxy with a significant amount of rotation (i.e., a disc). In addition, 
processes that feed star formation through the streaming of low-entropy 
gas along filaments (cold accretion) will also introduce rotation. Rotating 
structures resembling discs are found one way or the other over a wide 
range of cosmic distances, although disc galaxies at high redshift appear 
dynamically hotter than disc galaxies observed locally.” 
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Early-type galaxies 


Elliptical galaxies (E) and lenticulars (SO) constitute the other funda- 
mental group, termed ‘early-type galaxies’ (ETGs). These systems are 
mostly hot dynamical systems (even lenticular discs feature prominent, 
dynamically hot bulges). We will see below that their spheroidal mor- 
phology is supported by orbital anisotropy. Ellipticals are characterized 
by their apparent flattening (although note that the shape is dependent 
on projection), ranging from EO to E7, where the number corresponds to 
10 x (1—b/a) with b being the projected short axis and a the projected 
long axis. Overall, their stellar content is dominated by a roughly homo- 
geneously old and metal-rich population, suggesting an early and rapid 
process of star formation. This result led to the belief that massive ellipt- 
ical galaxies could form from a monolithic collapse of gas.* Lower mass 
ETGs have younger, metal-poor populations and may even feature faint 
levels of residual star formation — following the general trend of downsiz- 
ing, where the bulk of star formation progresses from massive to low-mass 
galaxies with cosmic time. 

The standard view of the formation of elliptical galaxies invokes a 
merging process where the progenitors can be either ellipticals or disc 
galaxies,* but should have comparable mass (major mergers). Additional 
growth mechanisms can be invoked via minor mergers, i.e., the infall of a 
low-mass satellite, but these processes are not expected to trigger a mor- 
phological change (although they can heat up the distribution of orbits, 
as in, e.g., the thick disc of the Milky Way). This merging scenario for 
the formation of ellipticals is consistent with the morphology-density rela- 
tion,” whereby elliptical galaxies are preferentially found in high density 
environments. The structure growth in these environments evolves at a 
faster rate in comparison with the field; therefore, the prevalence of ellip- 
ticals reflects the role of interactions among structures as a way to explain 
their formation. Compact massive galaxies are already found at redshifts 
z~ 2, i.e., when the Universe had roughly one-fourth of its present age, 
and are viewed as the progenitors of the cores of massive ellipticals and 
bulges at present time. However, the details of the formation process, in- 
cluding the role of major and minor mergers, and the initial cold accretion 
phase are poorly understood. 


5.2 Scaling relations 


The observable properties of galaxies have strong correlations that indic- 
ate, to a lowest order, a similar formation process driven by a reduced 
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Figure 5.1 Tully-Fisher relation from the SPARC disc galaxy sample. 
(Source: adapted from Lelli, McGaugh & Schombert, 2016, ApJ, 
816, L14.) 


number of parameters, most notably the mass of the galaxy (either total 
or stellar) or the average velocity dispersion of the stars. A blind search 
for correlations among galaxy properties via principal component ana- 
lysis found that many of the properties can be determined from a single 
parameter!° 


The Tully-Fisher relation 


One of the most important scaling relations in disc galaxies is the correl- 
ation between total luminosity and rotation velocity (figure 5.1).” The 
velocity that defines this relation is usually taken in the flat portion of 
the rotation curve, away from the centre. The scatter of this correla- 
tion is small enough to make it a standard candle and has been used 
to measure cosmological distances. The Tully-Fisher relation (TFR) can 
be determined from a simple dynamical argument invoking the virial 
theorem: 


M 
Vent = Sit E vior (5.2) 


However, this result implies that the TFR requires constant mass-to-light 
ratio across all disc galaxies, as well as constant surface brightness. The 
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Figure 5.2 Fundamental plane of a sample of early-type galaxies 
compiled from the Sloan Digital Sky Survey. (Source: adapted from 
La Barbera et al., 2010, MNRAS, 408, 1335.) Used by permission of 
Oxford University Press. 


latter is found to apply to the subset of high surface brightness galaxies 
(Freeman’s law). The relation depends on wavelength, increasing the 
power index of vyor towards the redder bands. Notice that at longer (i.e., 
redder) wavelengths, the light from the stellar component is mostly con- 
tributed by low-mass (mainly older) stars. Hence, variations in Y across 
the sample are much smaller than at bluer wavelengths. The correla- 
tion with respect to luminosity features a break at low-luminosities, with 
galaxies falling below the relation defined for the more massive galaxies 
(i.e., lower luminosity at fixed vpor). The correlation is extended to lower 
mass galaxies when using total baryonic (stellar plus gas) mass rather 
than luminosity, showing that the trend could, in principle, be fully re- 
lated to the baryonic mass; this result has been used as an argument to 
propose alternative models of gravitation that require no dark matter.® 


The Fundamental Plane and its projections 


The Fundamental Plane (FP) is a three-parameter correlation involving 
the size, velocity dispersion and surface brightness of early-type galaxies. 
These estimates are typically averaged inside an aperture extending over 
the effective radius or a fraction thereof. The data span a plane on this 
three-dimensional parameter space (see figure 5.2). By use of the virial 
theorem, it is possible to define a plane that appears to be tilted with re- 
spect to the observed FP. This mismatch can be related to several factors. 
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(1) A change in the mass-to-light ratio, Y « M® (see exercise 5.1), can ex- 
plain this tilt. Such an increase in M/L can be driven by a gradual increase 
of the contribution from dark matter towards more massive galaxies. 
(2) Rotation will “divert” part of the kinetic energy into bulk rotation, 
so that a systematic change in v/o with galaxy mass could introduce a 
similar effect. (3) The simplified model presented in exercise 5.1 also as- 
sumes that all galaxies feature the same surface brightness profile (i.e., 
their light distribution is homologous). A systematic trend, for instance, 
a change of the Sérsic index with galaxy mass, will produce a tilt of the 
plane defined by the simple model based on the virial theorem. The actual 
reason for the tilt of the FP is arguably a combination of the three.’ 


Exercise 5.1 


The Fundamental Plane is an important scaling relation of elliptical 
galaxies, where the average surface brightness (£), the size (R) and 
the velocity dispersion (o) are related by the equation 


Reo 


(these numbers have been slightly changed to make the problem 
consistent as is). Using the mass-to-light ratio Y =M/L to con- 
vert light into mass, show that the virial theorem leads to Rœ 
o? LIY, 

Let us now assume that the mass to light ratio changes system- 
atically from galaxy to galaxy, such that Y « M”. Find the value of 
a that brings the virial expectation in line with the observations. 


In addition to the three-parameter relation spanned by the Fundamental 
Plane, projections defined by just two parameters also provide interest- 
ing scaling relations of early-type galaxies, although not as tight, by 
definition, as the FP. 
1. Faber-Jackson relation: This trend involves the total luminosity 
(i.e., a combination of surface brightness and size) and the velocity 
dispersion: !° 

Lao, (5.3) 


and a typical value of the power law index is a~4. It is the equivalent 
of the Tully-Fisher relation in disc galaxies. Naively, one could guess this 
relation by assuming that while the dominant contribution to the kin- 
etic energy of discs comes from rotation, the kinetic energy in early-type 
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galaxies appears mostly in the form of velocity dispersion. However, there 
is a significant variation of the degree of rotational support in ETGs that 
correlates with mass, so this trend is not so trivial. 

2. Kormendy relation: This is the ‘easiest’ scaling relation from the ob- 
11 as it compares size (R) and surface brightness 
(x). Therefore it does not require any spectroscopic data. It can also be 


servational point of view, 


expressed by a simple formula: 
xR’, (5.4) 


and y is approximately —3, i.e., the surface brightness decreases with 
galaxy size. This scaling relation deeply encodes information about the 
dynamical history of ETGs. In simple words, the Kormendy relation states 
that bigger galaxies (broadly more massive galaxies) are ‘fluffier’. 


Exercise 5.2 


Using exercise 3.7, explain how a scenario that postulates the 
formation of early-type galaxies only via mergers can explain the 
Kormendy relation. 


The colour-magnitude relation 


In addition to the above relations, which depend mainly on galaxy dynam- 
ics, there is a scaling relation pertaining to the properties of the stellar 
populations. The colour-magnitude relation is especially significant in 
early-type galaxies, mostly lying on the red sequence (see figure 1.5). 
The more massive galaxies have redder colours. This relation has been 
exploited, for instance, to detect galaxy clusters, as they host a large frac- 
tion of elliptical galaxies. The redness of the colour of a population can 
be explained either by its (older) age, (higher) metallicity or (dustier) 
interstellar medium. The latter can be ruled out as elliptical galaxies do 
not have much gas (or dust) in them. A combination of the first two 
factors is thought to explain this trend, so that the star formation history 
of early-type galaxies is closely related to its total mass (more accurately 
its velocity dispersion). A detailed analysis based on spectroscopic obser- 
vations from the Sloan Digital Sky Survey revealed a strong correlation 
between the age!* and the metallicity!® of a galaxy with respect to its 
mass. This overly simple trend gets more complicated when dealing with 
colour gradients. In addition to morphology and the global dynamical 
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state of a galaxy, radial gradients provide further information about the 
assembly of galaxies, but this is something beyond the scope of this 
book. 


5.3 Rotation versus ‘pressure’ in early-type galaxies 


The tensor virial theorem, presented in section 3.9, allows us to relate the 
contribution from rotation and random motion to the shape of the galaxy, 
which, in turn, affects the gravitational potential energy. If the system is 
in equilibrium, there is no change in the moment of inertia tensor, and 
the virial theorem can be written 


QT jx + je + Wik = 0, (5.5) 


where the kinetic energy is split between a streaming motion component 
(T) and a random motion component (I). Using Cartesian coordinates 
and adopting the z direction as the axis of rotation, we can write 


Ta =Ty, 
Ty=0 if i Fj, 


(5.6) 


and likewise for TI and W. Hence, we only have two nontrivial equations 
for 5.5: 


2T yx + Tho + Wx = 0 
; (5.7) 


2Tzz T Izz T Wzz =0 


We can therefore relate the kinetic and potential energy terms as follows: 


9”, 6 


2Ty~ + Myx Wx 
2Tzz + Mzz | Wes 


where the final expression in this equation is derived by taking into 
account a spheroidal distribution of matter.'* Let us consider the two 
possible options that can explain the flattening of an elliptical galaxy: 
1. Rotation: If rotation dominates the kinetic energy budget and fully ex- 
plains the flattening of ellipticals, we can assume the velocity dispersion 
tensor is isotropic: 

Ilx = My = Mzz = Moĝ, (5.9) 


where oo is the mass-weighted velocity dispersion of the random motion, 
along the line of sight to the galaxy. The streaming component of the 
kinetic energy is 
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1 1 
Tx + Ty = 5 [ evaydx= 5M. (5.10) 


and vâ is the mass-weighted mean square rotation speed. Equation 5.8 
gives 


0 = Ar 3; (5.11) 
00 b 

2. Anisotropic velocity dispersion: The flattening can also be caused by 
a difference in the distribution of orbits (dynamical systems remember 
their past history, in contrast to a system in thermodynamic equilibrium). 
In this scenario, the streaming term vanishes: 


Tax = Tyy = Tzz = 0, (5.12) 


and we divide the random motion between a contribution on the rotation 
plane, and perpendicular to it: 


Thy = My =Mo2 
we ER (5.13) 
Izz = Mož 
giving 
b 0.45 
#-(?) : (5.14) 
Ox a 


In the extreme case of an E7 elliptical, 1 — b/a = 0.7, a full rotational sup- 
port would require v/o ~ 1.96, i.e., a rather large rotation velocity that is 
not observed in ellipticals, whereas a mild anisotropy (07/0; ~ 0.58) can 
explain the same shape. Figure 5.3 shows the relation between v/o and 
ellipticity in a sample of nearby early-type galaxies, where the solid dots 
represent slow rotators (i.e., hot dynamical systems) and the open sym- 
bols are fast rotators. The dashed line traces our simple expression from 
equation 5.11. A more accurate expression of the rotational support limit 
is shown as a solid line, roughly enveloping all the observed data points. !° 


5.4 A brief introduction to spiral arms in disc galaxies 


Disc galaxies often feature spiral arms. In fact, one of the main cri- 
teria for the classification of galaxies in the Hubble tuning fork diagram 
(see figure 1.3) is the morphology of these spiral features. From the 
observational point of view, the nature of spiral structure is correlated 
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Figure 5.3 Rotational support of early type galaxies as a function of 
ellipticity. (Source: data from the ATLAS* sample, Emsellem et al., 
2011, MNRAS, 444, 888.) Used by permission of Oxford University Press. 


with the global properties, such as the ratio of gas to stellar mass. Gas poor 
discs (lenticular galaxies: SO) have no spiral structure. Most spirals fea- 
ture two arms with conspicuous dust lanes on their inner edge. There are 
two main types of spiral morphology: grand design, with arms that can be 
traced over a large radial range, and flocculent, showing patchy, discon- 
tinuous spiral features (figure 5.4). The star formation rate is significantly 
higher in the spiral arms than in the rest of the galaxy. A consequence of 
this is the prominent display along the spiral arms of hot (OB, i.e., very 
young < 10 Myr) stars, as well as ionized gas (HII regions). 


The winding-up paradox 


A first guess at explaining spiral features would resort to differential rota- 
tion in the disc. We saw in chapter 4 that a significant part of the rotation 
curve of our Galaxy is flat. Hence the orbital speed is independent of 
radius, and the angular speed is Q (R) œ 1/R. This will imply that an exten- 
ded burst of star formation would gradually twist into a spiral feature as 
the central regions will wind faster than the outer ones. Consider a stripe 
of stars, all at the same azimuthal angle ¢9 =0 at time t=O. At a later 
time, the angular distribution will be #(R,t)=Q(R)t. The differential 
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Figure 5.4 Spiral arm morphology: flocculent (left) versus grand 
design (right). (Source: courtesy NASA/ESA HST and the Hubble 
Heritage Project.) 


Figure 5.5 Determining the pitch angle of a spiral feature. The dashed 
line follows a circumference, for which the pitch angle is i= 0. 


rotation will therefore shear the straight line into a spiral with a pitch 
angle (i), defined as the angle between the tangent to the spiral and the 
circle crossing at that point (figure 5.5): 


coti= rit =Rt 
=R p| 


dQ 
a —2At, (5.15) 


where we have made use of Oort’s A constant (equation 4.4). If the pitch 
angle is small, we can relate it to the separation, AR, between arms. Take 
two adjacent arms at the same azimuthal angle: 


dQ 
dR 


27nR = 27 
coti t 


I aR 
— At’ 


Qn = |2 (R+ AR) - Q(R)|t > AR = (5.16) 
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where the last step assumes AR «R. The separation decreases monoton- 
ically with time. Plugging in expected values for the Milky Way (t~10 Gyr; 
At = 151.36, R= 8.5 kpc), we find AR ~ 0.2 kpc, which is in conflict with 
the observational constraints. This result corresponds to a very small pitch 
angle (i~ 0.2°), whereas the observed value is around 5° for the earlier- 
type systems (Sa), increasing to 10° — 30° for the later-type ones (Sc, Sd). 
Hence, either the spiral arms are not long lived or the spiral pattern is not 
directly linked to differential rotation. 


Spiral pattern rotation 


A possible explanation for the presence of long-lived spiral arms invokes 
a resonance mechanism, such that stars are preferentially found near the 
spiral pattern, but, on average, the composition of stars within a spiral 
arm changes with time. An analogous example, closer to home, is the pile 
up of cars on a congested motorway. As cars approach the traffic jam, they 
will slow down, increasing the local density of cars. On the other side of 
the congestion, cars will accelerate, decreasing the density as they move 
away. If we look from above, we will see an overdensity of cars in a seg- 
ment of the motorway, but different cars join and leave this overdensity 
with time. Is it possible to produce a similar pattern in galaxies? From 
chapter 4, we know that stars in a disc galaxy do not move along circular 
orbits. They follow instead an elliptical orbit (an epicycle) that trails along 
a guiding centre representing the simple circular orbit (i.e., the location of 
the dLSR). The epicyclic motion has frequency « = 2,/—B(A — B), in con- 
trast with the angular frequency of the dLSR (Q). Let us assume there is a 
pattern of stars in a galaxy, i.e., an inhomogeneous distribution of stars 
(say a set of spiral arms), which moves with a constant pattern speed 
Qp. If the number m = (Q — Qp)/x« is a rational number, the trajectories 
will be closed and may allow for this density pattern to be sustained. The 
choice m= 2 is especially relevant, as it creates a two-armed structure. 
For a given choice of the pattern angular speed, there are three import- 
ant radial positions (figure 5.6, right panel). (1) At the position where 
Qp = (called corotation, CR), stars move with the pattern. Therefore, 
from a frame of reference rotating with angular speed Qp, we will see the 
stars simply tracing an epicycle. (2) At Qp = Q — «/2 (Inner Lindblad Res- 
onance, ILR), a star will move twice round its epicycle for every rotation 
around the galaxy, overtaking the pattern. (3) At Qp = Q + «/2 (Outer 
Lindblad Resonance, OLR), the star moves twice its epicycle for every 
rotation, falling behing the pattern. At any other radii, the orbits are un- 
closed trajectories. The left panel of figure 5.6 shows the radial variation 
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Figure 5.6 Illustration of Lindblad resonances. (Source: background 
image of M51 from the Digital Sky Survey.) 


of angular speed for an isothermal density profile; other profiles will look 
similar. Notice the Q — « /2 curve is flat outside of R ~ 4 kpc. If the pattern 
speed is similar to this value, the galaxy will extend its ILR over a large 
radial range; i.e., it would be possible to have oval-like orbits at many 
radii - when viewed from a frame of reference moving with the pattern 
speed. 

Figure 5.7 shows two different cases where an extended ILR can 
produce a pattern in solid body rotation. The figure nests orbits similar 
to the ones described above for the ILR over a range of radii. We expect 
these orbits to rotate with constant angular speed Qp. In the first case 
(left panel) the orientation is constant from the inside out, producing 
a bar-like feature that rotates rigidly. The second case (right panel) in- 
troduces a gradual rotation of the orientation of the orbits, producing 
a clear spiral feature, also expected to rotate rigidly. This explanation is 
purely kinematic, based on the motion of stars. It is meant only to motiv- 
ate the case for a resonance that can cause spiral features. Density wave 
theory details the mechanism by which the interaction of gas and stars 
in a disc galaxy can give rise to spiral arms in solid body rotation.'° The 
source of energy to drive the wave is still a matter of debate, and tidal in- 
teractions with a nearby object or the presence of a rotating bar may be 
involved. Note, for instance, that galaxy M51, shown in the right panel of 
figure 5.6, is tidally affected by a nearby galaxy (to the top of the figure), 
and many other grand-design spirals also have nearby galaxies that could 
be potential interlopers. We have considered in the example above a 
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Figure 5.7 Rigidly rotating features expected when the ILR is exten- 
ded. Either a bar (left) or a spiral (right) is produced, depending on the 
radial variation of the orientation of the orbits. 


two-armed spiral by choosing m = 2, but in general we could consider a 
pattern with m > 2. Stars and gas clouds will pass through the pattern 
with frequency m[Qp — Q(R)] perturbing the system. These motions with 
respect to the spiral pattern will affect the gravitational potential of the 
disc. It is found that stars respond so as to strengthen the spiral only if the 
perturbing frequency is slower than «. Hence, a spiral density wave can 
propagate only in the region where this condition holds, i.e., 


Q-—K/M SQ <QA+K/mM. (5.17) 


The extrema of this interval are the inner and outer Lindblad resonances. 
Patterns with many arms (m > 2) will therefore span a narrow region in 
the diagram (e.g., the interval between Q — «/4 and Q+«/4 will shrink 
in figure 5.6), impeding the formation of extended spiral structure with 
many arms. The vast majority of grand-design spiral galaxies display two 
arms. Note that at the corotation radius, stars move along with the pat- 
tern, whereas inside (outside) corotation they will move faster (slower) 
than the pattern. When gas clouds hit the spiral pattern, they get com- 
pressed, triggering star formation. This explains why most of the young 
stars are found along the spiral arms. Once they form, the stars over- 
take the spiral feature (inside corotation). Hence, we should expect an 
offset between the gas, ongoing star formation, and more evolved stars, 
following the (orthogonal) distance to a spiral arm. This is one of the 
characteristic properties exploited to test density wave theory.!” Gas is 
more affected by the spiral pattern than stars, because of its smaller 
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random motion. As stars move in the galaxy, they increase their random 
motion. The distribution of older stars in the disc is more homogeneous 
(i.e., less affected by the presence of spiral arms). Therefore, lenticular 
galaxies (SO, disc systems with very little gas) have no spiral structure, 
whereas later-type systems (Sb, Sc) normally feature the most prominent 
arms. 

Nevertheless, at present, numerical simulations cannot create a 
spiral pattern with solid body rotation, as predicted by density wave the- 
ory. Recent research suggests that spiral arms are short-lived, in corota- 
tion with the stars. '° This scenario would be compatible with the presence 
of flocculent arms, traditionally explained by the onset of detonation 
waves caused by star formation, propagating around the disc and get- 
ting sheared by differential rotation. However, grand-design spirals pose 
a challenge to numerical models when tidal interactions are not invoked. 
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6 
Galactic chemical enrichment 


Most of the chemical elements in the periodic table were synthesized in 
the cores of stars. The Big Bang, although very hot during a brief period 
of time, could not keep the required condition of temperature and dens- 
ity for long enough, resulting only ina 25 per cent production by mass 
of helium from hydrogen, along with traces of a few heavier nuclei. The 
synthesis in the Universe of elements such as C, O, Mg and Fe can take 
place only in the central regions of stars, where the conditions for thermo- 
nuclear reactions are adequate. Therefore, the build-up of the elements 
in galaxies — observed in the stellar and gaseous components - reflects 
the past star formation histories inside galaxies. Tracing the composition 
of gas and stars provides a powerful tool for understanding galaxy form- 
ation. This chapter presents an overview of chemical enrichment along 
with the differential equations that govern the variation of elemental 
abundances. The simplification of the instantaneous recycling approxim- 
ation allows us to produce very simple solutions that give insight about 
the different scenarios of galaxy evolution. We relate these solutions to 
the observed scaling relation between mass and metallicity and the solu- 
tion of the G-dwarf problem in Milky Way stars. Abundance ratios are also 
introduced as a powerful discriminant of star formation histories. 


6.1 Nucleosynthesis and the formation of galaxies 


The formation and evolution of galaxies can be split into two complement- 
ary channels. One deals with the growth of the dark matter halos within 
which galaxies live. This growth — the mass assembly history (MAH) - is 
driven mainly by gravitational interactions, from the growth of the small 
density fluctuations at early times (chapter 7) to merging processes that 
gradually create more massive structures. The second channel involves 
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the star formation history (SFH) of the galaxies inhabiting these halos 
and can be studied via galactic chemical enrichment. As star formation 
proceeds, the chemical elements synthesized in stars are progressively 
incorporated into subsequent stellar generations, although a fraction of 
these elements, along with the H+He gas, can be expelled from galaxies. 
A simple description of these proceses allows us to understand the mech- 
anisms underlying the transformation of gas into stars. This is arguably 
one of the most complicated problems in astrophysics, as a large number 
of physical processes play an important role, including the hydrodynamics 
of gas inflows and outflows, the cooling and fragmentation of gas clouds, 
feedback from star formation and activity from an Active Galactic Nucleus 
(AGN). 

Central to the analysis of galactic chemical enrichment (hereafter 
GCE) is the concept of stellar yields, defined as the mass of a given 
chemical element synthesized in a star of mass m. The yields are cal- 
culated from models of stellar structure and evolution, aspects that are 
beyond the scope of this textbook. However, the most fundamental as- 
pect is the significantly different role played by massive stars (>8Mo). 
They contribute the majority of chemical elements, and feature compar- 
atively short lifetimes (25 Myr for a 10 Mo star, in contrast with 10 Gyr 
for a 1 Mo star). Therefore, chemical enrichment closely traces the star 
formation history. The presence of an additional nucleosynthetic channel 
involving a binary system — type Ia supernovae — introduces an additional 
timescale that will be presented at the end of this chapter. The specific de- 
tails of a star formation region can, in principle, leave an imprint on the 
detailed chemical composition of the stars formed within it. The so-called 
chemical tagging analysis is a promising method by which GCE could help 
identify the dynamical origin of a population of stars in our Galaxy (or 
nearby resolved systems). 


6.2 General aspects of galactic chemical enrichment 


In the simplest form, the equations of GCE trace the evolution of the gas 
and stellar mass, as well as the amount of chemical elements. Often, we 
add together all chemical elements — except for H and He - defining metal- 
licity (Z) as the mass fraction in these elements with respect to the total. 
As a reference, our Sun has a metallicity Zo ~ 0.02. An alternative estim- 
ator of metallicity uses an individual element, typically iron, with [Fe/H] 
defined as the logarithmic ratio, by number, of iron atoms, with respect to 
the solar value. In general, one can find stars and stellar populations with 
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very low metallicity ([Fe/H] <—3, and even lower), but high metallicity 
populations are not found above a factor of ~2 — 3 times the solar value. 


The initial mass function 


Of all parameters describing a star, its mass is the most important one 
(Russel-Vogt theorem). Given a stellar mass, we can infer its lifetime, 
yields, evolutionary phases and remnant properties. Additional factors 
such as chemical composition and rotation will also affect these, but we 
assume that, to a good approximation, only mass defines lifetime and 
yields, essential parts of the GCE equations. Therefore, a fundamental 
aspect of chemical enrichment is the distribution of stellar masses in a 
population at birth. This is the initial mass function (IMF). Local ob- 
servations of stellar populations support the idea that this function may 
be universal, although recent evidence in star-forming and quiescent 
galaxies reveals potential variations, possibly towards a top-heavy IMF 
in star-forming systems (i.e., an excess of massive stars with respect to 
the standard IMF); and a bottom-heavy IMF in massive quiescent galaxies 
(i.e., an excess of low-mass stars with respect to the standard). Neverthe- 
less, we will assume here a unique definition of the IMF, given by a single 
power law: 

o(m)= ay aml, (6.1) 

dm 

where T = 2.35 defines the Salpeter IMF. We also need to consider the 
mass range of the IMF, i.e., the interval of mass within which stars can 
be found. The low-mass end is caused by the threshold in the onset of 
thermonuclear reactions in the stellar core, at around 0.08Mg. Lower 
masses do not have enough gravitational energy to achieve the high cent- 
ral temperatures needed. Instabilities in the stellar atmospheres control 
the upper limit, although the actual value is not clearly defined. A value of 
100 Mo is commonly adopted for this limit. Quite often, the IMF is defined 
on a logarithmic scale in mass: 


hee cence (6.2) 
dlogm 


with the Salpeter! IMF corresponding toa slope u =T — 1 = 1.35. We also 
need to use the normalization mass scale for a given IMF, defined as 


Mhigh 
my = f m@(m)dm. (6.3) 


Mow 
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This normalization is typically fixed at a solar mass, my =Mo, and some 
textbooks redefine the normalization, removing the mass units from this 
expression by choosing my = 1. Note that in this case, the IMF is defined 
per unit mass, instead of per number of stars. There are a number of addi- 
tional functional definitions of the IMF, where the main effect is either to 
taper off the low-mass end with either an additional power law (or sets of 
power laws), as in the Scalo or Kroupa definitions, or to replace the power 
law at low masses with a lognormal distribution, as in the Chabrier IMF. 


Exercise 6.1 


Show that for a Salpeter IMF, defined in the mass interval 0.1 
-100 Mo a supernova is expected per 184 Mo of stars formed (as- 
suming that only stars with mass m > 10Mo undergo a supernova 
explosion). The star formation rate of our Milky Way galaxy is 
approximately one supernova per century, and the average star 
formation rate is 1 Mọ yr_!. Which way should we change the IMF 
to reconcile these estimates? 


Remnant mass and returned fraction 


Once a star reaches its endpoint (after a time tm), we assume that a rem- 
nant is left behind: a white dwarf for low- and intermediate-mass stars 
(m <8-10M,) and a neutron star or black hole for massive stars. This can 
be quantified by defining a remnant mass (wm) that remains locked away 
at later times. The rest of the mass is ejected back into the interstellar me- 
dium (ISM), including both ‘gas’ and ‘metals’, which become available for 
the next episodes of star formation. Therefore, if we consider the form- 
ation of a single population of stars (all created at the same time), we 
can define the fraction in mass returned back to the ISM. This so-called 
returned fraction can be related to the IMF as follows: 


1 Mhigh 
R(t) = — J m) (m — Wm)dm, (6.4) 


MN me 


where m; is the stellar mass corresponding to a lifetime t. This function 
is zero at t = 0, gradually increasing with time towards an asymptotic 
value at late times, when stellar lifetimes are as long as the age of the 
Universe. Typical values of the returned fraction at late times are ~0.3 for 
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the Salpeter IMF and ~0.5 for the Chabrier IMF. Since the latter features 
a lower fraction of low-mass stars — which effectively lock up most of their 
mass forever — the amount of gas available for recycling is significantly 
higher. 


Exercise 6.2 


Compute the returned fraction at late times of the Salpeter IMF (ad- 
opt 1 Mo as the low-mass threshold of the integral), if the remnant 
mass of a star with initial mass m is given by (in solar mass units): 


wm  [0-L(m/Mo)+0.5, m<10Mo, 
Mo 


1.5,  m>10Mo. 


Stellar yields 


Another important component of GCE equations is the yield (pm), defined 
as the mass fraction in a given element, that is synthesized by the star 
(of mass m) and returned to the interstellar medium once it reaches its 
endpoint, after a time tm. An equivalent version of the returned fraction 
is the net yield, defined as 


1 Mhigh 
Yp= TER | mpm (m)dm. (6.5) 
=R Jm 


Star formation rate 


The star formation rate (y) is an essential function in GCE equations, as 
it drives the transformation from gas into stars. The Schmidt law assumes 
that y is some power law of the gas mass available for star formation. A 
simple argument gives 


poa BË > yO =ke). (6.6) 


tdyn 
From an observational point of view, only projected surface mass 
densities can be measured. The equivalent version of the above on two- 
dimensional surface formation rates and surface gas mass densities is the 
Kennicutt law: 


GALACTIC CHEMICAL ENRICHMENT 


111 


112 


Dy (t) = kgr Di O. (6.7) 


Exercise 6.3 


Assuming LxM?-> — consider only the main sequence and neglect 
evolved phases of stellar evolution — compare the average mass and 
average luminosity of a population at birth and after 1 Gyr (the life- 
time of a 2.5 Mo star is 1 Gyr). Neglect the contribution to the mass 
from remnants. 


6.3 Basic equations of galactic chemical enrichment 


The basic equations of GCE? describe the evolution of the mass in stars 
(Ms) and gas (Mg), as well as the amount of “metals” in the gas phase, 
described by the metallicity (Z). We assume the total mass is M=M, + 
Mg. More advanced equations separate the gas component into several 
‘phases’, mainly a cold phase (which controls the star formation rate) and 
ahot phase. Furthermore, the metallicity can be split into the contribution 
from individual elements, or from groups closely related to their synthesis 
reactions, such as the a elements, iron-peak elements, or s- and r-process 
elements. Gas flows are simplified by two time-dependent functions, an 
infall rate, f(t) and an outflow rate o(t). The basic equations describing 
the evolution of the mass components are 


dM 

af? 

dM, _ 

ae =yp-E : (6.8) 


dM, 

g 
-Aan E+f—o 

dt YE+f 
The quantity E(t) is the ejection rate, i.e., the amount of gas contributed 
by stars at the end of their lives. Noting that wm is the remnant mass, we 
can write 
Mmax De 
E(t) =| (m—Wm) W(t — tm) (m)dm, (6.9) 
mt 

where yy = y/my is the star formation rate per unit mass, so defined as 
to give E(t) the units of a mass per unit time. The evolution of the metal 
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content is given by 


d(ZM, 

o = Fh Epi =o: (6.10) 
The enrichment process is described by a homogeneous mixture of gas 
and metals, with instantaneous mixing. The ejection of metals is given by 


EzO= [0m —Wn)ZCe— m) + mpm] Öl- tmm, (6.11) 


mt 


where pm are the stellar yields, i.e., the mass fraction of a star of mass m 
converted into metals, and ejected at the end of the stellar life. 


The instantaneous recycling approximation 


These equations can be simplified if we assume that stars are divided 
into two classes: those that live forever (masses below some threshold 
m < mo) and those that die out as soon as they are born (m > mg). This is 
called the instantaneous recycling approximation (IRA). The approxima- 
tion is well justified because most of the yields are produced by M = 8Mo 
stars, which have much shorter lifetimes than ~Mo stars, which lock most 
of the mass into stars. The ejecta (see equations 6.9 and 6.11) are now 


E(t) =Ry (0), 
Ez(t) =RZOY O +y -R [1-20] vO. 


(6.12) 


The factor [1 — Z(t)] in the second equation is ~1, given that the metal- 
licity is always Z « 1. The equations of chemical enrichment in the IRA 
become 


dM, _ 

q Ry 

Te = AR +f- o . (6.13) 
ae) =-Z(1-R)W+y(1—R)w +Z;f — Zo 


Instead of tracing the mass in metals (ZM,), one can write the equation 
for the evolution of the metallicity (Z), as follows: 


dZ 
q YAW + G-a. (6.14) 


Msg 
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Neglecting gas infall and outflow, we see that the total amount of metals 
ever formed at some time t is 


t Mmax 
ZM: + ZM = | | mpm — m$ mdm | ds =yMs. (6.15) 
0 


Ms 


The closed box model 


The above equations can be applied to basic models of chemical evolu- 
tion to derive analytic expressions of the stellar metallicity distribution. 
The simplest option is the closed box model, which has no infall or out- 
flows, with initial conditions M,(t = 0) = Mo; Z(t = 0) = 0. The total mass 
remains constant M=M;+M,=Mp at all times. Integrating the IRA 
equations gives 

Mg Mg 


dZ dM 
= —Z=ylnu;!, h =—-_=—. (6.16 
3; M; ylnyg where jg Mg +M; Mo ( ) 


Note this expression is valid as long as Z < 1, so it will break down at late 
times, as ug — 0. The closed box model produces a significant amount of 
low-metallicity stars, which is as expected since the star formation rate 
at early times (i.e., with low-metallicity gas) is very high. Comparisons of 
this model with observed stellar metallicities in the solar neighbourhood 
show a large mismatch, termed the ‘G-dwarf problem’. We can quantify 
this difference by estimating the mass fraction in stars with metallicity 
below some threshold Z. We need to rewrite the IRA equations, relating 
the metallicity to the stellar mass: 


dZ dM; 
y E Mo — Ms 


= M,(<Z)=Mo (1-70). (6.17) 


In order to derive this expression, we use the fact that the metallicity in the 
gas phase increases monotonically in this model, so that the stellar mass 
when the gas phase metallicity is Z comprises only stars with metallicities 
below this value. We can quantify, for instance, the mass fraction in stars 
with metallicity below one quarter of the solar value. Firstly, we need to 
derive the yield, by use of the observational constraint of the gas mass 
fraction in the solar neighbourhood, at present time: 


E 5Mo pc~? 
gas oo OE 6a (6.18) 


NOW) = = 
Hyl } Xgas+stars IOMo pc~? 
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and the metallicity in the gas phase is at present Z ~ Zo. Therefore, 
from equation 6.16 we find y = 0.43Z 6, which, when plugged into equa- 
tion 6.17, gives a mass fraction Ms(< Zo/4) ~ 0.45. The observational 
constraint is ~0.02, i.e., a substantially smaller fraction. Therefore, the 
closed-box model does not give a valid representation of the solar neigh- 
bourhood as it overproduces low-metallicity stars. 


The open box models: Infall and outflows 


To mitigate the G-dwarf problem posed by the closed box model, we 
need a mechanism to reduce the production of low-metallicity stars. The 
simplest solution beyond the closed box model is to consider that stars 
are being formed from the infall of gas supplied by an external reservoir. 
By assuming a small initial gas mass, we prevent the early locking of a 
high-mass content in stars when the metallicity is too low. To simplify the 
model, let us assume no outflows, and infall perfectly balancing the star 
formation rate, such that the gas mass is kept constant: Mg = Mo => f = 
(1—R)y, leading to 


M,(<Z)=Mo n( 2), (6.19) 
y-Z 


where we also assume that the gas reservoir has zero metallicity. As in the 
closed box case, the net yield can be derived from 


dZ dM; Z 


1 
=a =< s eliz, (6.20) 


Therefore, in the solar neighbourhood (ug=0.1) we get y~Zo, and 
M;(< Zo/4)=0.032, in line with the observational constraints (see 
figure 6.1). 


Exercise 6.4 


Show that the closed box and infall models presented above lead 
to the following distribution of stars with respect to metallicity, as 
shown in figure 6.1: 


dN te $, closed box, 


x 
dloge |, infall, 
1-¢ 


where ¢=Z/y. 
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Figure 6.1 Comparison of the observed metallicity of nearby stars 
from the Geneva-Copenhagen Survey", shown in grey, with the closed 
box and infall model presented here (the latter labelled ‘Infall #1’). A 
second, more realistic distribution from a numerical model of chem- 
ical enrichment? is shown as ‘Infall #2’. Used by permission of Oxford 
University Press. 


Holmberg et al., 2007, A& A, 475, 519. 
>Ferreras, Wyse & Silk, 2003, MNRAS, 345, 1381. 


It is also interesting to incorporate the effect of outflows in the chem- 
ical enrichment process. For instance, let us assume that the outflow 
rate is proportional to the star formation rate, defining o=cM,, where 
c is the proportionality constant. This ansatz is justified by the fact that 
supernovae-driven winds, triggered effectively by the instantaneous star 
formation rate (SFR) when IRA is adopted, drive the gas outflow rate. If 
there is no infalling gas, and the initial condition is Mg =Mo, and zero 
stellar mass, we get 


TOE i hes -1 
Z=7 7 he =Yer IN ug ; (6.21) 


equivalent to the closed box model, where the effective yield is the net 
yield, reduced by a factor (1 +c). Likewise, the cumulative stellar mass 
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below some metallicity threshold is 
M a7, 
M;(<Z)= 7 [z-e va], (6.22) 


again, equivalent to the closed box model. By increasing c, one could bring 
this mass fraction in line with the observations. 


Exercise 6.5 


Consider now the leaky box model, with the same properties and 
initial conditions as the closed box model, but with additional infall 
(f) and outflow (0) terms, where f =o=cM, (i.e., the infall and 
outflow rates are balanced and correspond to a fraction, c < 1, of 
the rate of change in stellar mass). Derive the mass fraction in low- 
metallicity (Z < Zo /4) stars, and contrast with the observed ~0.02 
fraction. 


The outflow model can also be presented in the light of the observed 
mass-metallicity relation (figure 6.2), where the most massive galaxies 
feature the highest metallicities. If we make a simple relation between 
the observed gas-phase metallicity and the effective yield, a functional de- 
pendence would be established between c —i.e., the gas outflow efficiency 
-and galaxy mass. Such a relation could be motivated by an escape velo- 
city argument, where the gravitational potential well would ‘modulate’ 
the amount of gas ejected in outflows. 


Exercise 6.6 


A galaxy fuels its star formation via infall, but without any out- 
flowing material. At time t=0 the gas mass is Mo and the stellar 
mass is zero. The infall rate is at all times a fixed fraction (aw) of the 
star formation rate (f =ay), and the star formation rate follows 
a simple law: y =kMg. Find the time evolution of the gas frac- 
tion u = M,/(Mz+M;). What is the asymptotic behaviour of u as 
t + œ? What happens if a = 1? 
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Figure 6.2 Observed correlation between the stellar mass of the 
galaxy and the (gas-phase) metallicity, measured with respect to the 
oxygen abundance. In these units, solar metallicity corresponds to a 
value of 8.7. (Source: adapted from Tremonti et al., 2003, ApJ, 613, 
898.) 


6.4 Chemistry as a cosmic clock 


So far, the GCE equations adopt a single variable to describe metallicity 
(Z). The next step would entail separating the production of elements into 
a elements and Fe-peak elements. The main reason for this choice is that 
the former are mainly produced by core-collapse (Type II) supernovae, 
which are the final evolutionary stage of massive (28 Mo) stars. In these 
stars, most of the iron-rich core is trapped into a remnant (neutron star or 
black hole), locking this material away. The Fe yield of Type II supernovae 
is generally rather low (<0.1 Mo per supernova). In contrast, Type Ia su- 
pernovae - triggered in close binary systems involving at least one white 
dwarf -release a significantly higher amount of Fe (~0.7 Ma). Since such 
a scenario requires the presence of a white dwarf, a delay is expected in 
the production of Type Ia supernovae with respect to Type II, and thus 
iron is incorporated later in subsequent generations, as long as star form- 
ation is present. Although estimating this delay is a complicated task that 
depends on a large number of factors involving the production and evol- 
ution of binaries and the accretion of material in a way that would yield 
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Figure 6.3 Correlation between the fraction of old stars (vertical axis, 
given as the mass fraction in stars older than 10 Gyr) and the chemical 
abundance ratio [#/Fe] (where 0 corresponds to the solar value). Note 
that higher ratios are associated with galaxies where the old component 
is more prominent, suggesting a relation between the star formation his- 
tory and chemical enrichment. (Source: adapted from de la Rosa et al., 
2011, MNRAS, 418, L74.) Used by permission of Oxford University 
Press. 


a supernova explosion, we can approximate this delay by a Atra, which 
allows us to use the abundance ratio [a/Fe] as a cosmic clock: galaxies 
with high stellar [w/Fe] result from an intense, short-lived episode of star 
formation that shut off quickly, before Type Ia supernovae had any chance 
to contribute to chemical enrichment. If [w/Fe] is low, the star formation 
rate was more extended. For most models of the Type Ia progenitor, the 
delay time (Atra) ranges between 0.5 and 2 Gyr. For example, massive 
early-type galaxies feature metal-rich and old stellar populations, with 
substantially higher [a /Fe] than stars in the solar neighbourhood, reflect- 
ing a very different star formation history (figure 6.3). The old and high 
metal content suggests an early and efficient formation process to quickly 
build up the gas-phase metallicity, and to avoid locking too many low- 
mass (i.e., long-lived) stars at low metallicity. The high [a/Fe] reinforces 
the need for a high formation efficiency, requiring most of the stars to 
be formed within <1-2 Gyr. Massive galaxies are observed at high red- 
shift, potentially the cores of high-mass early-type galaxies at present. 
They have a stellar mass of about 10!!Mo at redshift z~2. These systems 
must sustain a very intense formation scenario, where the average star 
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formation rate needs to be about two orders of magnitude higher than in 
the Milky Way at present. In the Milky Way there is also a signature of 
nonsolar abundance ratios, as low-metallicity stars (formed in the early 
stages of evolution) have high [w/Fe], a property used, for instance, to 
disentangle the populations of thin and thick disc stars (see section 4.1). 


Notes 


1 Salpeter, 1955, ApJ, 121, 161. 
2 Galactic chemical enrichment was a field pioneered by Beatrice Tinsley (Tinsley, 
1980, Fundam. Cosmic Phys., 5, 287). 
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/ 
The growth of density fluctuations 


Galaxies start as minute density fluctuations in an otherwise homogen- 
ous and expanding Universe. During the first phases of galaxy growth, 
the fluctuations were so small that simple equations can be solved to 
trace this early evolution (linear phase). An overview of cosmology is 
necessary, including the treatment of the expanding Universe through 
Friedmann’s equations, although a simple Newtonian argument will suf- 
fice here, in lieu of a more complex treatment with general relativity 
(outside the scope of this textbook). The concept of critical density and 
the evolution of the background matter and radiation densities will be 
presented, leading to the equations of an expanding fluid with a small 
perturbation, parameterized with the density contrast. The Einstein—de 
Sitter case will be presented as an easy-to-solve model that gives an ac- 
curate representation of the Universe during the main phase of linear 
growth. A comparison with the density fluctuations at present in the 
Cosmic Microwave Background leads to one of the most robust proofs 
that dark matter must be present in large quantities in the Universe. The 
linear phase is followed by a more complex stage of nonlinear growth 
that can be simplified with a model of spherical collapse and virializa- 
tion, also presented here. Once the fluctuation is in virial equilibrium 
(forming a so-called dark matter halo), we consider the issue of gas 
cooling, responsible for the typical sizes and masses of galaxies. The 
distribution of the dark matter halos leads us to a statistical treatment 
based on the Gaussian distribution, following the Press-Schechter ar- 
gument. This treatment results in a hierarchical buildup of structure. 
Finally, galaxy clustering is presented as a way of exploring the underlying 
cosmology. 
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7.1 A cosmology primer 


Galaxy formation, in its earliest stages, involves the growth of small dens- 
ity fluctuations in an otherwise homogeneous and isotropic Universe. This 
approximation, called the cosmological principle, represents the central 
tenet of our current paradigm of cosmology, allowing us to use a highly 
simplified metric as the solution to Einstein’s equations of general relativ- 
ity. This textbook avoids a relativistic treatment, as stellar and galactic 
velocities are, at most, ~10° kms~! «c, and the ratio of gravitational 
to rest-mass energy is W/Mc? ~ 10~°M10/Rkpes where Mọ is the mass in 
units of 10!°Mo, and Rxpe is the typical size of the stellar system in kpc. 
However, it is necessary to invoke the spacetime metric to quantify the 
(important) role of the expansion of the Universe in the growth of density 
fluctuations. The Friedmann-Lemaitre-Robertson-Walker metric (FLRW) 
depends on only two parameters, an overall dimensionless scale factor, 
a(t), which accounts for the physical separation between two galaxies 
with cosmic time, and the curvature constant, «, with three different 
values: —1, 0, +1 for a closed, flat and open Universe, respectively: 


d 2. 
ds? = cdt? — a?(t) | ” 4r (a9? + sin? aag?) | . (7.1) 
1— «r2 


The scale factor is set to unity at present time: a(to) = 1. Observations of 
the angular distribution of temperature in the Cosmic Microwave Back- 
ground (CMB) impose stringent constraints on the curvature, leading to 
a flat Universe: « =0. Although we cannot measure the scale factor dir- 
ectly, it is possible to relate it to a direct observable: the wavelength of 
a photon emitted from a distant source (say at cosmic time t). Since all 
scales are affected by a(t), we can write 


A(t) _ Alto) 


— E i 7.2 
a(t) (to) ve 


where A(t) = Ao is the rest-frame wavelength, and A(to) = obs is the ob- 
served wavelength. Noting that the redshift (z) is defined as: 1+ z= 
Aobs/A0, we find that 


a(t) = (7.3) 


1+z° 


Therefore, a galaxy at redshift z=1 lives in a Universe whose size is 
half of the size of the present one. We will see that in cosmology one 
can interchangeably use time, the scale factor or redshift to track cosmic 
evolution. 
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Another important parameter is the rate of change of the scale 
factor. It can be measured as the recession velocity between distant galax- 
ies. Note that in order for these estimates to be valid within the context of 
the cosmological principle, we need to work with large enough distances 
so that homogeneity and isotropy hold. Therefore, the relative velocity 
between two nearby galaxies will be affected by local interactions (termed 
‘peculiar velocities’). Incidentally, regardless of the expansion of the Uni- 
verse, our nearest neighbour, the Andromeda galaxy, is moving towards 
us! Itis for much more distant galaxies that the effect of peculiar velocities 
is negligible, and the cosmological expansion (Hubble flow) dominates. 
A fundamental relation is the Hubble law, where the recession velocity 
r) is 

Vp = U9 4 = H(Od, (7.4) 

a(t) 

and d is the separation between the galaxies. Since all galaxies are affected 
by the cosmic expansion, it is sometimes useful to define a comoving dis- 
tance, recom = d/a(t), thereby effectively factoring out the expansion.! In 
comoving coordinates, the only motion between two galaxies is caused 
by peculiar velocities, a useful framework when tracing the growth of 
density fluctuations. The Planck 2015 value (see table 7.1) of Hubble’s 
constant is H(to) = (67.8 + 0.9) kms~! Mpc~!. Noting that peculiar ve- 
locities among galaxies can be as high as 10° kms~!, we can infer that 
distances between galaxies much larger than ~20 Mpc are needed to be 
able to ‘feel’ the cosmological expansion. 


Friedmann’s equations 


Friedmann’s equations provide a solution to Einstein’s General Relativistic 
(GR) equations in the FLRW metric. Although GR is beyond the scope of 
this textbook, we need to consider GR to effectively link the geometry of 
spacetime (the metric) to the distribution of matter and energy in the sys- 
tem. Friedmann’s equations allow us to relate the evolution of the scale 
factor, a(t), to the (homogeneous) density in its various guises: radiation, 
matter or dark energy. We will show here a simplified derivation based 
on Newtonian mechanics which arrives at the correct equations. Let us 
model the Universe as a sphere with radius r, expanding with the Hubble 
flow. The force on a test particle with mass m is 


GM(<r)m 4r Gm 
r2 = 


mř = rp(t) äs- paw, (7.5) 


THE GROWTH OF DENSITY FLUCTUATIONS 


123 


124 


where we have identified the radius of this ‘sample Universe’ as the scale 
factor. Since the Universe is an isolated system, conservation of mass 
gives pa? = poa} = constant. The subindex zero refers to the present time. 
Therefore, multiplying equation 7.5 by å gives 


åä = 5 ae ac (poa) å. (7.6) 
Integrating this equation leads to 
a= T pa? — Ke, (7.7) 


with the last term being an integration constant. In GR this term has the 
meaning of the constant curvature allowed by the FLRW metric. To de- 
rive the correct equation from this Newtonian approximation, we need to 
include the contribution to the energy from pressure. Following a simple 
thermodynamic argument, consider the amount of heat flowing into the 
Universe, with volume V. Noting that the internal energy is € = Vpc*, we 
have 


AQ= AE + pAV=VA(pc2) + (pc? +p) AV. (7.8) 


Since the Universe is a closed system, AQ=0, and since the volume is 
3 


Vaca’, we can write AV/V = 3Aa/a. Therefore: 
a (pd) 0 (7.9) 
dt `a cay i 


Taking the time derivative of equation 7.7, and using 7.9 to describe 
dp/dt, we arrive at the second of Friedmann’s equations, namely: 


@=-al[p +35]. (7.10) 


However, this Newtonian example cannot describe the so-called cos- 
mological constant (A). Friedmann’s equations including this constant 
are 


82G 1 
& = = pa? kc + 3 Aa? 
(7.11) 
is 4nG 3p 1 
= A 
a 3 a [v + Z | + 3 a 
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Dimensionless densities 


The first of Friedmann’s equations can be written in terms of the Hubble 
parameter: 


8rG kÊ A 
H?(t) = — . 12 
(t) gi ats (7.12) 


This equation introduces a density scale, the critical density, that is time 
dependent: 
pc(t) = 3H? (t) /8xG. (7.13) 


At present time, the critical density is p¢(tg) = 1.879 x 10-7°h? kg m-3, 
or, in more convenient units, 2.775 x 10!4h#Mo Mpc~°, where Hubble’s 
constant is usually defined as H(to)=Ho=100hkms~!Mpc™! (so, 
roughly, h=0.7). The density is usually split between radiation (,) 
and matter (pm). These two components vary with a as py xat and 
Pm x a°, respectively. The —3 exponent is a volume dilution factor as the 
Universe expands (from conservation of energy) whereas the —4 in the 
radiation field takes into account the additional loss of energy because 
of the redshift (note that the energy of a photon is e, =hv =hc/A). We 
also define the equivalent, dimensionless density parameters, taking the 
(time-dependent) critical density as reference: 


t 
wo- O 
za 7.14 
Pmt) ‘ ) 
Qm(t) = 
Perit) 


The cosmological constant can also be associated to an additional energy 
component (dark energy), with density parameter: 


Qa (7.15) 


3H5 
Plugging these values into the Friedmann equations at present time 
(a(to) = 1) allows us to write the curvature as 


tmt- -1 


CHA j (7.16) 


which means the Universe is flat (zero curvature) if the sum of the dens- 
ity of all components is the critical value. A simple case often used is the 
Einstein-de Sitter model, where 9, = 24 = 0 and 9m = 1, which can be 
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Figure 7.1 Evolution of the density parameter of the radiation (Q,), 
matter (Qm) and cosmological constant (Q,), as a function of the scale 
factor. The grey shaded area marks the region 2m > 0.9, where the Uni- 
verse is well described by an Einstein—de Sitter model. (Cosmological 
parameters taken from table 7.1.) 


shown to have a simple solution for the scale factor of 


t 2/3 


The evolution of the density parameter is shown in figure 7.1, which 
illustrates the regime where an Einstein—de Sitter model gives a valid rep- 
resentation (grey shaded region). Note this interval corresponds to the 
major phase of galaxy formation. 


Exercise 7.1 


Show that the dimensionless matter density scales with redshift as 


Qmo(1 +2)? 


Qm(Z) = 3 
aE E ETE ETEN, 


where the ‘0’ subindices refer to the values at present time, and 
Q, is the corresponding term for the curvature. In our standard 
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cosmological model (ACDM), the cosmological parameters are, 
approximately Qm,o ~ 0.3; Qy,0 * 1074; 24,0% 0.7, so that the 
curvature is zero. Show that during an earlier epoch, the Universe 
was closely represented by an Einstein—-de Sitter model. 


The age of the Universe 


One can estimate the age of the Universe using the following integral: 


t | "y f worn dt fe f o (7.18) 
= — a= = P ` 
T Jo o da o aHa) Jo +z)H@) 


Note we can use the time coordinate t, the scale factor a or the redshift z to 


track ‘cosmic time’. Theorists usually prefer the scale factor, whereas ob- 
servational cosmologists prefer redshift. Applying Friedmann’s equation 
to an Einstein—-de Sitter model gives 


H(z) =Ho(1+2)*””, (7.19) 
leading to 
1% dz 2 
t= = = 52h 1G ~9.3 Gyr. 2 
” n (1+z2)5/2  3Ho 6 yr ~ 9.3 Gyr (7.20) 


The most general case can be written as the equivalent of equation 7.18, 
using 7.12 and 7.13 to write Hubble’s parameter: 


H(2) = HE"? (2) = Hoy mo +2? + 2y 01 +2)4 + 2,01 +2)2+ 24,0, 
(7.21) 


and so: 


pes f -n (7.22) 
Ho Jo CFDEZ®O 


which, for a vanilla-flavoured ACDM cosmology, with Qm,o = 0.3, Qa,0 = 
0.7, 2,,0 ~ 0, and, therefore, zero curvature, gives an age of the Universe 
at present of: ty = 9.43 h-t Gyr~ 13.7 Gyr, i.e., significantly older than 
the Einstein—-de Sitter case. To determine the age at an arbitrary redshift, 
z, we simply have to replace the lower limit of the integral in equation 7.18 
by z. 
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Distance(s) in cosmology 


Distance is a fundamental parameter in the interpretation of the observa- 
tional properties of galaxies. In this succinct primer, we need to refer only 
to a fundamental difference between distance indicators. The distance 
directly derived from the metric (equation 7.1) is the comoving radial 
distance. If we consider the radial trajectory of a photon from a distant 
galaxy, at coordinate r, and our observing point, at r = 0, the distance is 


4 dr a=1 ceda Z g 
D.= See = dz, (7.23) 
: i; VIr? [ a?H(a) I H(z) 


and the expression for Hubble’s parameter with respect to redshift is 
given by equation 7.21. Note that for cosmologically small distances, 
we can approximate this equation by cz ~ HoDe, i.e., Hubble’s law of 
expansion. 

However, this distance is not what we would use when measuring 
the angular extent of a galaxy. In this case, we measure an angle (0) given 
by the ratio of its (comoving) size and the comoving distance to us. If the 
physical diameter of the galaxy is ø, the comoving size is 6/a=9(1 + z); 
therefore 


= Ø +2) d , where Da(z)= aal 
D(z) Da(z) 1+z 


(7.24) 


and we define the angular diameter distance Da as the one that relates the 
physical size with the angular extent of the galaxy. Moreover, we also ob- 
serve the flux from galaxies as the luminosity (L) per unit area. A sphere 
centred at the position of the observed galaxy, touching our observing po- 
sition, has area 42 D2. However, the photon rate will decrease by a factor 
a~! = (1 + z) because of the cosmological time dilation, and the energy of 
the photon will decrease by an additional factor of (1 + z) because of the 
cosmological redshift. Therefore, the observed flux can be written 


T 
7 4x D2(z) T 4nD2(z)(1 + 2z)2 
Lo 
= 4nD2(@)’ where Dj(z) =D;(z)(1+2), (7.25) 


where Lo is the intrinsic (i.e., rest-frame) luminosity of the galaxy, and 
D, is defined as the luminosity distance. Note that all three distances 
converge at low redshift, where large-scale effects become negligible. 
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Exercise 7.2 


Show that while the surface brightness of galaxies (section 1.2) 
remains constant when the distance is small, over cosmological dis- 
tances it varies as £ (z) œ (1 + z)~4. This is called ‘Tolman dimming’. 


7.2 Linear regime 


A perfectly homogeneous Universe can be easily described by the equa- 
tions presented in the previous section. However, such a density distribu- 
tion would not develop structure. Hence, we need to introduce minute, 
but nonnegligible perturbations in the density field at early times. An 
inflationary scenario posits that microscopic fluctuations (over Planck 
scales, i.e., driven by the quantum nature of gravity) will expand at early 
times in an exponential way, moving outside the horizon — i.e., the re- 
gion within which perturbations can be in causal contact. Being outside 
the horizon means the different areas of the perturbation will not be in 
causal contact. Therefore, these fluctuations are frozen (they cannot in- 
teract), and require general relativity (outside the scope of this book) for 
a correct description. Later on, because of the steady expansion, the ho- 
rizon catches up with the extent of these fluctuations. For the sake of a 
simple, pedagogical argument, it is assumed that a standard Newtonian 
treatment in an expanding background suffices. 

Let us consider the equations that describe the evolution of a density 
fluctuation as a fluid. First of all, we write down the unperturbed solution, 
(Vo, po, Po, Po), where we adopt the Lagrangian notation (i.e., the time 
derivatives follow the fluid) d/dt = 3/3t + ©- v): 


d 3 
o =—p0V -Vo 
dv 1s z 
W app = Vint (7.26) 
dt po 
Vo = 4r Gpo 


These three expressions describe the conservation of mass, the force (i.e., 
Euler) equation and Poisson’s equation, respectively. These equations are 
now modified with small perturbations in the velocity, density, pressure 
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and gravitational potential fields: (ôv, 5p, dp, 6®). Retaining the terms to 
the lowest order, we get 


d (de Ba AR: 
aan, (pees Yee, ve 
a) ” 


HBV) o eS is 7 
+ (8- V)¥p = ——V (8p) — V(S) 
po 


dt 


(7.27) 


VŽZ(8 D) = 4r Gop 


The next step involves a change of coordinates from physical distance 
(x) to comoving distance (r), factoring out the effect of the cosmological 
expansion: 

x=a(t)r 


d 
je 2h ka i ean} A (7.28) 
dt dt 


p(X, t) = pot) [1+ 67,0] 


> 


Physical velocities (v) are decomposed into a Hubble flow (Hx) and a 
peculiar velocity (ù). The density is also split between a smooth, back- 
ground average (po) and a density contrast (5). Note that the comoving 
separation between two galaxies that move only with the Hubble flow 
(i.e., ‘at rest’ in an expanding Universe) remains constant with cosmic 
time. Finally, we arrive at the linear perturbation equation for the density 
contrast: 


$ 2 
+2 ($) i- 4rGm= $v +00,3. (7.29) 
a a 


This equation has an oscillator term (« 5) driven by the background mat- 
ter density, a damping term caused by the expansion (« 6), and a driving 
force, on the RHS, with a characteristic scale given by the sound speed 
(cs = ./6p/5e). This equation is similar to the growth of a fluctuation in 
standard Newtonian dynamics (say, the collapse of a large cloud of gas 
in a galaxy), with the difference being the damping term caused by the 
expansion. This term fundamentally affects the rate of structure growth. 
For example, in an Einstein-de Sitter Universe (Qm = 1), and assuming 
dust-like matter (c; = 0), 


3 23 à 2 
a(t) = (3tot) see 
2 a 


87G 
re eel AU 
3H2 


(7.30) 


2 
=1 > 4rGpo = =5, 
ET 
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giving the following linear perturbation equation: 
3 s 2 
ô+ —ô— -ô =O. (7.31) 


If we try a solution of the type ô(t)œ t”, we find solutions for n= —1 
(decaying mode) and n = 2/3 (growing mode). Neglecting the decaying 
mode, we find that 


2/3 


b(t) axt”? «a(t. (7.32) 


Therefore, a density fluctuation will grow linearly with respect to the scale 
factor, in contrast with the density perturbation in a gas cloud within a 
galaxy, which follows an exponential rate: the cosmological expansion 
dampens this growth. 


Exercise 7.3 


Neglect the pressure term in the linear growth of density fluctu- 
ations (i.e., the term with V2), and consider an empty Universe 
(i.e., zero matter and energy density). This is the Milne model. 
Show that in this case, density fluctuations will remain constant 
with time. 


Jeans mass 


The treatment of the density growth presented above traces only the 
time evolution, but it does not take into account the scale of the fluc- 
tuation. The linear density growth equation (7.29) can be applied to 
a specific lengthscale. One can describe the spatial distribution of the 
density fluctuation as a Fourier series in space, 


b= >) i OCN (7.33) 
k 


where the coefficients, 6, quantify the contribution of lengthscale à = 
27/k to the fluctuation. For simplicity, we will assume spherical sym- 
metry, considering k a scalar quantity. In the linear regime, each com- 
ponent, 6,, evolves independently, giving 


be +2 (5) bp = (47Gpo — cok?) Sk (7.34) 
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where we have used the comoving wave-number k, = k/a. We can solve 
this equation with the ansatz 


ôk, (t) = 5k.,0€ (7.35) 
which leads to a dispersion relation 
w =c2k? — 4r Gpo. (7.36) 


Hence, the small density fluctuations will oscillate if the RHS is positive. 
If it is negative, the fluctuations will have exponential growing/decaying 
behaviour. This means density fluctuations will only grow if 


27 x \i/2 
G 2k? > A> dy = =e, | — : 3 
4r Gpo > cs ké = A> Az = «(Z) (7.37) 


The second part refers to the size of a (comoving) fluctuation (å), defining 
the Jeans length. Any density perturbation whose size is greater than A; 
will grow. This threshold can also be given as a mass (called the ‘Jeans 


mass’): 
4r 


M;= 3 


por3. (7.38) 
The Jeans instability reflects a balance between pressure forces (related 
to the sound speed) and gravitational forces (from po). If a fluctuation 
involves a mass greater than Mj, then gravity will overcome pressure, 
and the fluctuation will collapse. For smaller masses, sound waves will 
counteract the gravitational pull, preventing the formation of a structure. 
Let us consider the growth of galaxies after photon decoupling. At 
this epoch (z < 1,100), we can make the simplifying assumption that the 
sound speed of the gas is given by that of monoatomic hydrogen gas: 


| SkpT 
Cs = 2B a37? kms7!, (7.39) 
3mp 


with T3=T/1000K. Close to decoupling, T3~+3, and the comoving 
Jeans length is 0.014 Mpc, leading to a Jeans mass around 1.6 x 10°Mo 
(where we made the assumption of an Einstein-de Sitter Universe and 
Ho = 70 km s Mpc~!). Therefore, anything more massive than a globular 
cluster can, in principle, form a collapsed structure. However, the coup- 
ling of the ordinary matter to the photons at earlier times washed out any 
potential fluctuations over these scales. We will see in the next section 
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how an additional matter component is needed to produce the high levels 
of density contrast we see in the Universe today. 


The need for (nonbaryonic) dark matter 


The linear behaviour of a perturbation implies that a long time is needed 
to evolve from the the early linear regime (6 >> 1) to the present distribu- 
tion of matter in the Universe. Furthermore, observations of the Cosmic 
Microwave Background — which give a glimpse of the density fluctuations 
in the baryonic content of the Universe at photon-electron decoupling, 
i.e., 380,000 yr after the Big Bang — reveal temperature inhomogeneit- 
ies around T/T ~ 107°. Considering adiabatic temperature fluctuations, 
i.e., keeping the entropy, S œ ag ig constant, we can write 


ôT 1ên, 1p 


Tam 36 (7.40) 


It is assumed that before recombination, photons and electrons/baryons 
were well coupled (via Thomson scattering). Hence, the measurements of 
the CMB temperature fluctuations imply a density contrast at decoupling 
(at redshift zcmg ~ 1, 000) around ô(tcmg) ~ 3 x 107°. Using a simple ar- 
gument for the growth in an Einstein-de Sitter model (equation 7.32), we 
find those fluctuations should be at present time as follows: 


a(to) 
a(tcms) 


5(to)~ (teu) ( ) =at + cmp) ~ 0.1. (7.41) 


This means the contrast of the baryon density field at the time of de- 
coupling is not large enough to create galaxies today! Hence, we need 
an additional component, not coupled to photons, so that the growth of 
their density seeds took place at earlier stages, unimpeded. After decoup- 
ling, the baryons fall on to the (larger) density fluctuations of this (dark) 
matter field, allowing for larger values of the density contrast. 
Therefore, let us consider a two-fluid component involving baryons 
and dark matter (the contribution from radiation is negligible at decoup- 
ling and afterwards). The growth of density fluctuations of both fluids 
is described by the linear perturbation equation (neglecting the pressure 
term): i 
ög + 2 bp = 41 G(p0,B5B + /0,DM5DM) 
i (7.42) 
dpm + 2Žġpm = 4r G(p0,pM5pm + (0,B5B) 
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Figure 7.2 Evolution of the density contrast with cosmic time, 
parameterized with the scale factor, a. (Plot produced with the 
PyCAMB code of Lewis & Challinor, Astrophysics Source Code Library, 
ascl:1102.026.) 


Assuming a dominant contribution from the dark matter, we find (for the 
Einstein—de Sitter case) that 


dpm = dpm(O)a(t), (7.43) 
which goes back into the first equation, giving 
bp +2 ($) g z= 4r Gapo,pmôpm (0). (7.44) 


Rewriting the differential equation with respect to the scale factor (or 
redshift), instead of time, we find that 


1 
Baim € = =), (7.45) 


where zy is some reference value (obtained from the constant of in- 
tegration). Hence, at high redshift, during the earliest phases of evolu- 
tion, z —> zy, and the baryon fluctuations are negligible (e.g., during the 
epoch of recombination/decoupling, zy ~ 1,000), whereas later on, at 
z X 1, 000, the density fluctuations of the baryons follow the dark matter 
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distribution, ôg ~ édpy. Figure 7.2 shows the evolution of the density 
contrast of the dark matter (5,), baryons (6,) and radiation (6,,) as a func- 
tion of the scale factor. A standard ACDM model is assumed, with a scale 
k = 0.3 Mpc—!. Note that all fields evolve similarly during the radiation- 
dominated phase (a ô «a? trend is shown as reference). At later times, 
the photon/baryon mix undergoes a set of oscillations, and at decoupling 
(a~ 107°), the baryon matter follows the evolution of the cold dark mat- 
ter density fluctuations. At later times, the trend 5 «a is also shown as 
reference. 


7.3 Spherical collapse 


The spherical collapse model is the simplest way to follow the evolution of 
a density fluctuation from its initial growth in the linear regime (6 « 1) - 
where it ‘detaches’ from the Hubble flow - to the final state of a virial- 
ized, stable structure (i.e., a dark matter halo). We begin with a spherical 
region, with radius R(t), evolving within an expanding background. The 
force equation drives the evolution of the region: 


G 
oo a oo(1+8)R. (7.46) 


Notice the analogy of this equation with the evolution of the scale factor 
in an expanding Universe: 


= pa. (7.47) 


Originally, the fluctuation is so small that the region expands with the 
background, but the slight overdensity gradually slows it down. At some 
time, tmax, the sphere will halt its expansion and start collapsing (turn- 
around). Conservation of energy gives 


1 /dAR\? GM GM 
R R 


z \ ae — =€(M)= A (7.48) 


The time evolution is therefore 


1 1 \J77 


This equation can be readily solved if we define R = Ryax cos? «æ. Finally, if 
we redefine the angular parameter a =*/2 — ° /2, so that R =0 at 6 = 0 and 


THE GROWTH OF DENSITY FLUCTUATIONS 


135 


% 
Q, 
“3 
Be 
Zoe 
i) 
4 
E 
% 
A 
i TEE i! uiul 


100 : >= >= a 

= 10 =E -z 

Si E 3 

SS E J 

3 L Turnaround - 

2 TF os 

E / a 

O1b 4 
0.01 


TM 


0.2 0.4 0.6 
Scale factor a 


Figure 7.3 Schematic figure depicting the evolution of the density ina 
fluctuation following the spherical collapse model. 


R=Rwmax at 0 ="/2, we arrive at the following parametric equations: 


R(@) =Ro(1 — cos 0) 


t(0) = to(0 — sin 8) , (7.50) 
9(6 — sin 0)? 


2(1 — cos 0)? 


p (0) = po (0) 


which describe the trajectory of a cycloid. Note that the density equation 
takes into account the evolution of the background density in an Einstein- 
de Sitter model, i.e., poxa? and acxt?/3. In these equations, to = 
/ RRiax/8GM «./1/Gp, i.e., the dynamical timescale (see section 1.3). 
Therefore, denser fluctuations — on average found at earlier times — will 
have shorter collapse timescales. 

In this model (see figure 7.3), the fluctuation starts expanding with 
the background. However it gradually slows down until it stops expand- 
ing at the turnaround time, when 6 = x. At this time, the radius of the 
fluctuation is Rmax = 2Ro. The density at turnaround is 5.559. The spher- 
ical collapse model allows us to determine the epoch at which a structure 
collapses (virializes): 
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1+2max _ a(2tmax) 
l1+zvr  a(tmax) 


— 22/3 ~1.59, (7.51) 


(for an EdS model). Hence, if turnaround occurs at zmax = 10, the system 
is virialized at zyjp ~ 6. We can also probe the properties of the virialized 
halo by applying energy conservation: 


3GM? W 1 /3GM2 
E(tmax) = E(tvr) > “R ( ) 


5Rmax 2 2 (5RvR 
(7.52) 
1 
= RvR = JEMAX: 
leading to a halo density at collapse of 
PvR = 8p (tmax) ~ 44.4p0 (tmax) = 
(7.53) 


1+ ZMAX 
1 +ZvR 


3 
= 44.4 ( ) po(tvir) = 180po(tvir), 


providing a relationship between the density of a halo and its cosmic 
epoch of collapse: 


pvr ~ 180 perit 02m,0 (1 + zvr)”. (7.54) 


For instance, for the properties of the Milky Way (Mpm ~ 3 x 1011Mọ; 
Rv ~ 50 kpc) we find a virialization epoch of zyr ~ 3, whereas the col- 
lapse of a galaxy cluster (Mpm ~ 101°Mọo; Ry ~ 1.5 Mpc) happens at a 
later epoch zyr ~ 1. This result reflects the hierarchical nature of struc- 
ture formation, with small, denser fluctuations occuring first. However, 
this result concerns the formation of dark matter halos. The galaxies 
are actually formed through the subsequent collapse of baryonic mater- 
ial (H, He gas) into the centres of these halos. Unfortunately, a direct 
mapping of dark matter formation and galaxy formation is not straightfor- 
ward, and requires a proper understanding of processes involving, among 
others, the collapse of the infalling gas and its cooling, star formation and 
potential feedback effects. These complex mechanisms are loosely termed 
‘baryon physics’. Moreover the formalism assumes no mass mergers after 
the collapse of the halo. Nevertheless, it provides powerful insight about 
the relation between the epoch of formation and the properties of the dark 
matter halos. 
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Exercise 7.4 


A dark matter halo forms in an Einstein—de Sitter universe, col- 
lapsing at redshift z,i,, with a mass density profile given by 


2 
Ve 


Sy. r<ro 
4nGr2’ i 


plr)= 


0, r>Tro; 


where ro represents the radial extent of the halo, and ve is the velo- 
city of any circular orbit in the halo. Show that the radial extent of 
the halo is 

Ve 


ro = 32 
10H (Zyir) 


where H(Zyir) is the Hubble parameter at virialization. Assuming 
that the stars in a galaxy have the same circular velocity as its dark 
matter halo, how could we compare the formation epochs of two 
galaxies with the same ve? 


The linear density equivalent 


In the next section we will extend spherical collapse to a statistical de- 
rivation of the halo mass function. First, however, it will be useful to 
derive the equivalent density contrast for collapse in the linear regime, 
i.e., assuming that the difference between the density of the clump and 
the background density can be expressed to the lowest order. In this case 
it is useful to rewrite equation 7.50 with respect to the turnaround radius 
(Rmax) and time at turnaround (thax) as 


(7.55) 
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Writing out s as a function of t by retaining the lowest order terms in @ 
leads to 


4 0 (7.56) 


Note that the lowest order term gives the evolution of the size of the 
perturbation if it is moving with the Hubble flow (i.e., following an 
Einstein—de Sitter model). Now we write the density contrast in the lin- 
ear approximation as 1 + djin = Shack /S, Where the numerator corresponds 
to this lowest order behaviour that follows the background expansion. 
Retaining again the lowest order term leads to 


on (621)2/3 E co] 


3 
lin = géo = 1.69, (7.57) 


where the numerical solution corresponds to the linear prediction at 
collapse time, i.e., 7 = 2. 


Exercise 7.5 


Numerical simulations of structure formation suggest that dark 
matter halos follow a universal profile (Navarro-Frenk-White, or 


NFW): 
olx= r _ 50 Perit 
Ts x(1+x)2’ 


(assume spherical symmetry), where rs is the scale radius of the 
halo, perit is the critical background density, and 69 is the dens- 
ity contrast. The concentration of the halo is defined as c=rvypr/rs. 
Following the spherical collapse model, show that 


180 c? 
3 In(di+c)—c/(1+c) 


o= 


7.4 Press-Schechter formalism 


The combination of the linear density growth equation (7.29) and 
spherical collapse (7.54) allows us to trace the evolution of single density 
perturbations from the small seeds with ô « 1 all the way to the formation 
ofa stable halo. The next step involves dealing with the whole distribution 
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of density fluctuations as a random field, identifying the formation of 
halos in a statistical way. Such a treatment will enable us to define a mass 
function of virialized halos, much in the same way that we defined the 
luminosity function of galaxies in chapter 1. However, this process gets 
complicated as structures can form within structures: a galaxy in a cluster 
can be considered part of the dark matter halo that formed the progenitor 
galaxy early on, or part of the cluster that incorporated the galaxy/halo 
system at later times (cloud-in-cloud problem). We need a way to associ- 
ate a density fluctuation that eventually forms a halo with a (mass or size) 
scale. To do that, the density field is convolved with a spherical top-hat fil- 
ter, defining at each point the smoothed average density within a radius 
Ry as 


ô; Œ; Ry) = i ECOWE-— T; Rdr, (7.58) 


where the function W (f; Ry) is constant inside |r| < Rm, and zero outside. 
A given scale Ry corresponds to a mass M = 4r (p)R3,/ 3. The Press- 
Schechter formalism assumes that the probability that ês > ê: —i.e., the 
density contrast for collapse — is the same as the fraction of structures 
present at time t, contained in halos with mass greater than M. The linear 
theory equivalent of the density of a collapsed structure shown in equa- 
tion 7.54 corresponds to 5, = 1.69 (equation 7.57). The density contrast, 
6(r), can be described by a random Gaussian field, and likewise for the 
smoothed version, so that the probability of having a structure of mass 
greater than M is 


[0,6] 

—82 /20? (M) 
n e°s dô; x F(> M), (7.59) 
o(M)/2z Js, ceri 


where the variance of the Gaussian distribution is 


P[> del = 


o? (M) = (62(7;Ru)), (7.60) 


and F(> M) is the fraction of halos with mass greater than M. A problem in 
this formalism is that as M — 0, the mass variance ø? (R) diverges, and the 
total probability P — 1/2, as if only half of the total mass would eventu- 
ally collapse into halos. Therefore, only the mass contained in overdense 
regions, i.e., above average, are assumed to collapse into halos. However, 
underdense regions can form part of collapsed structures, and Press and 
Schechter argued,” without proof, that all matter will eventually collapse 
into halos, therefore requiring an extra factor of 2 to account for this: 
F(>M)=2P[> ê], so that the number density of halos with masses in 
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the range M and M + dM is 


n(M)dM = A 22) 8c 82/202 | dino 
M 0M am M2 o 


InM dM. (7.61) 


The Press-Schechter formalism gives a useful tool to map the evolution 
of structure in a statistical way. Only halos with mass M form in a sig- 
nificant number if o (M) > ô.. The time evolution of the mass function is 
implicit in the variation of o (M) with cosmic time, defined as ø (M, t) = 
o (M, 0)D(t), where D(t) is the growth factor of fluctuations, for instance 
in the Einstein—de Sitter case, D(t) = (t/tg)?/*. An alternative way to in- 
troduce the time evolution is to adopt a different density contrast for 
collapse with cosmic time, following 6,(t) =6,(0)/D(t). In the standard 
ACDM framework, o (M) decreases with mass, implying that halos form 
in a bottom-up way, from small, dense structures at early times towards 
massive halos at later times. If we assume a simple scaling relation for the 


variance, 
M —a 
= = .62 
we obtain the following halo mass function: 
M a—2 (M 2a 
n(M) =no a) e (if) j (7.63) 


where nọ is the normalization factor, and 


25 
M,=Mo| 5 (7.64) 
(z) 
is the characteristic mass scale. Note this result is reminiscent of the 
Schechter function (see equation 1.5), presented in the description of the 
luminosity function as a power law with an exponential cutoff. As cos- 
mic time increases, the density contrast threshold decreases as 1/D(t), 
pushing the exponential cutoff towards higher masses, as expected in the 
progressive growth of structure. 


7.5 Correlation function 
In addition to the census of halos with respect to mass, given by the 


Press-Schechter formalism, one can also consider the correlation of the 
density fluctuations: Are they randomly distributed in space, or do they 
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preferentially cluster? The clustering information provides a strong dis- 
criminant among different structure formation scenarios. We quantify 
this property with the two point auto-correlation function &(r), defined 
as the excess probability of finding a pair of galaxies (or clumps, or halos) 
separated by a radius vector r, with respect to a purely homogeneous dis- 
tribution. The number of galaxy pairs within infinitesimal volumes dV, 
and dV2 is therefore 


dNp(F) = NAIL + &(F)]dVidV2, (7.65) 


where N? gives the value for a homogeneous, uncorrelated distribution. 
This can be related to the density contrast: 


Np (7) = p(X)dVi p(X +7)dV2 = pll + EOL +647) ]dVi dVo. 


(7.66) 
Averaging over a large number of volume elements, 
Np (F) = poll + (8(%)5& +7))1dVidV2. (7.67) 
Hence, the two-point correlation function is 
E(r) = (88+ T). (7.68) 


Invoking isotropy, this function is assumed to depend only on distance but 
not on direction. The observed correlation function is found to be well 
represented by a simple power law, such that 


=y 
E(x (=) (7.69) 
ro 


over physical scales between 100 h~!kpc and 10h~!Mpc, with a charac- 
teristic length scale ro =5h~!Mpc, and y = 1.8. This trend depends on 
galaxy morphology, with elliptical galaxies being more clustered than 
disc galaxies. The connection between the correlation function and the 
power spectrum can be established if we take the Fourier transform of 
the density contrast: 


> V —ik t3 
O= ys | He : rq k, 


1 E (7.70) 
ô; = 7 [oër 
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Using Parseval’s theorem relating Fourier transform pairs, we find that 


1 2/2) 73, V 1273 
sf dPr= Ony3 ic. d'k, (7.71) 


where the mean square of the density fluctuations (LHS= (52)) is given 
with respect to the quantity 5,17, defined as the power spectrum P(k). 
For each value of k, the power spectrum gives the ‘intensity’ of dens- 
ity fluctuations over lengthscales à = 27/k. Using the definition of the 
autocorrelation function from equation 7.68, we get 


EO) f P(kekF ak. (7.72) 


(278 
Hence, the two-point autocorrelation function and the power spectrum 
are Fourier pairs. Assuming spherical symmetry, and inverting the previ- 
ous integral, we can compute the power spectrum from observations of 
galaxy clustering: 


SOUS dk (7.73) 
kr 


1 CO 
P(k)= = 
(k) v f &(r) 
The primordial power spectrum is assumed to be a simple power law: 
Po(k) xk"; (7.74) 


this implies a correlation function (for the primordial structures) of 


i k L i i 
gory | Peak ~r (R43) y M7033, (7.75) 


where we have transformed the distance r into a mass scale M via M ~ 
pr?. We can relate this scaling to the density fluctuations (eq. 7.68) over 
amass M: 


50(M) = (82) 1/2 x M7 @+3)/6 (7.76) 


At very early times, fluctuations grow as ô ~ a? (t) (see figure 7.2), and the 
amount of mass within the horizon goes as M x a? > a x M'/3. Hence, the 
density contrast of fluctuations over a mass scale M at early times — when 
they cross the horizon — is 


8(M, t) x a2(t)89(M) x M2/3M~+3)/6 x yr" /6, (7.77) 
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Figure 7.4 Transfer function for three different cosmological models: 
a standard ACDM scenario (see table 7.1); a baryon-dominated model 
where the Q, and Qm values of the ACDM model are swapped; and a 
hot dark matter (HDM) model, dominated by massless neutrinos. (Plot 
produced with the PYCAMB code of Lewis & Challinor, Astrophysics 
Source Code Library, ascl:1102.026.) 


Hence, for n=1, the fluctuations are independent of mass. This is a 
desirable behaviour that defines the Harrison-Zeldovich spectrum. The 
primordial density fluctuations will be affected by a number of physical 
processes, removing power from specific lengthscales. For instance, hot 
dark matter removes power from small scales, because of free streaming. 
The transfer function T(k) relates the original fluctuations to the evolved 
ones, namely: 


ôk ©) = de (O)T(K)D(O), (7.78) 


where D(t) is the growth of a fluctuation with cosmic epoch, derived from 
the geometry of spacetime, e.g., for an EdS universe D(t) = a(t) = t?/? (see 
equation 7.32). Some of the effects from the transfer function can be seen 
directly in the two-point autocorrelation function. The wiggles at scales 
~ 100h7! Mpc (related to the evolution of the growth function with time, 
as shown in figure 7.4) reflect the Baryon Acoustic Oscillations. They are 
avery powerful tool to constrain cosmological models. 


Types of dark matter and the power spectrum 


The fundamental aspect that defines the role of dark matter in the form- 
ation and evolution of structure is the average energy per particle at 
decoupling, i.e., when the expansion rate exceeds the interaction rate of 
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dark matter particles with the rest of the components, much in the same 
way as decoupling between the photon and electron/proton components, 
but at earlier times. We can differentiate three types of dark matter that 
will leave different imprints on the density distribution: 

Hot dark matter (HDM) is made of particles that are relativistic at de- 
coupling. They are low-mass (but non-zero), fast moving particles (e.g., 
standard neutrinos) that do not cluster into small structures, because a 
large mass (> 1014Mọ) is needed to keep them gravitationally bound. 
Since the speed of sound is close to the speed of light, the Jeans mass 
must be very high. Structure formation evolves from the top-down, i.e., 
starting with the formation over large scales/masses, with fragmentation 
later producing galaxy-type structures. 

Cold dark matter (CDM) consists of particles that are nonrelativistic at 
decoupling. They clump over all scales, leading to a bottom-up process 
of structure formation. Examples are heavy neutrinos, axions (these ones 
may be nonthermal relics), or the more generic WIMPS (Weakly Interact- 
ing Massive Particles), such as the neutralino (involving supersymmetry). 
Warm dark matter (WDM) has intermediate properties between HDM 
and CDM. The formation of small structures is suppressed, perhaps solv- 
ing an excess of small halos in a CDM scenario. However, this also delays 
the formation process of galaxies, which is at odds with the observa- 
tion of galaxies with very old stellar populations at present time, and the 
presence of massive galaxies at relatively high redshift. 

As the growth of density fluctuations as a function of mass de- 
pends on the nature of dark matter, the measurement of the density 
fluctuation power spectrum at a given cosmic epoch will constrain the 
nature of dark matter. The scale-dependent growth of fluctuations is de- 
scribed by the Transfer Function, shown in figure 7.4 for three different 
cosmological models. Note that the HDM model will suppress density 
fluctuations over small scales, with respect to the standard ACDM frame- 
work. The baryon-dominated scenario also suppresses fluctuations over 
small scales, because of its early coupling with the photon fluid (indeed, 
justifying the need for nonbaryonic dark matter; see section 7.2), and 
introduces a significant amount of oscillations. 

CDM is the favoured model — along with a cosmological constant, 
defining the ACDM model. The Planck 2015 best-fit cosmological para- 
meters (see table 7.1) show that only ~ 1/6 of the matter content in the 
Universe is in the form of ‘baryons’. You should also contrast this with the 
much lower fraction (~3%) expected in galaxies (figure 1.6), suggesting 
that a large amount of baryons was ejected from halos and should live in 
the intergalactic medium, in the form of a diffuse gas. Figure 7.1 shows 
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Table 7.1 Cosmological parameters. 


Parameter Symbol Value 

Matter density m 0.3089 + 0.0062 
Baryon density 2p 0.0486 + 0.0010 
Dark energy density Qa 0.6911 + 0.0062 
Hubble constant Ho 67.74 + 0.46 
Power spectrum normalizarion og 0.8159 + 0.0086 
Power spectrum index n 0.9667 + 0.0040 


Source: Planck collaboration, 2016, A&A, 594, 13. 


Note: All parameters are dimensionless, except for the Hubble constant, given in 
kms~! Mpc7!. 


that the expansion at present is dominated by the cosmological constant 
(dark energy). 


7.6 Cooling and the masses of galaxies 


Right after the formation (virialization) of a dark matter halo, we can 
assume that the gas is homogeneously distributed in this halo, at the 
virial temperature. However, the gas will lose energy, settling into the 
central region, cooling down and eventually forming stars. A very simple 
argument, invoking the cooling function, leads to the physical proper- 
ties of structures where galaxies form. The cooling function, A(T, Z), is 
defined such that the energy loss per unit time of gas with (electron) 
number density n is dE/dt = n* A(T, Z). Hence, one can define a cooling 
timescale as 3 nk 

tcool = g = oe (7.79) 
The cooling function is governed by processes that result in emission of 
photons, with the subsequent loss of energy from the system. This com- 
plex derivation is simplified by assuming a balance between the collisional 
ionization and recombination of the different atom and ion species. When 
a galaxy forms, two fundamental timescales are confronted: the cooling 
timescale and the dynamical timescale, where the latter depends only on 
the mass density. The collapse time of a homogeneous sphere with mass 
M and radius R is 


n is 


~5x 107 
i yr (G cm-3 


1 
tayn * = 7.80 
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Figure 7.5 Region where galaxy formation is allowed, following 
the Rees-Ostrikercriterion.* The solid and dashed lines delimit the 
parameter region corresponding to fixed values of mass and size, 
respectively. 


Rees & Ostriker 1977, MNRAS, 179, 541. 


The difference between dynamical, cooling (and cosmological) time- 
scales determines whether a baryon gas can form a galaxy. There are 
three main regimes to consider: 

1. Hp | < teol? No structure can form, and the gas stays as a diffuse 
component. 

2. Ho 1S teool > tayn: A hot structure can be formed, but the cooling is not 
efficient enough, forming a hot, pressure-supported gaseous halo, as in 
the intracluster medium (see chapter 9). 

3. Ho 1> tayn > teool: Cooling proceeds, triggering the collapse of the gas 
in the central region of the halo, leading to galaxy formation. 

Therefore, the parameter that controls galaxy formation is R= 
teool/tayn- Only for R <1 can the gas cool and form stars. This can be 
directly related to the typical sizes and masses of galaxies. Figure 7.5 
shows the galaxy formation regime within the grey shaded region — inside 
which R < 1, representing gas with solar chemical composition. The black 
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wiggly line — corresponding to zero metallicity — would delimit instead a 
narrower region for star formation, as expected, since a higher metalli- 
city leads to a more efficient cooling process. The dashed and dotted lines 
mark the trends of fixed mass and radius, respectively, assuming an ideal 
gas, and a 10 per cent contribution in baryons. 


Exercise 7.6 


Show that the slope of the lines in the log-log plot of figure 7.5 for 
structures with constant radius have slope 1, and the lines at con- 
stant mass have slope 3. Adopt an ideal gas equation of state and 
impose hydrostatic equilibrium to find the answer. 


Notes 


1 Note that the FLRW metric (equation 7.1) is defined with comoving coordinates. 
2 Press and Schechter, 1974, ApJ, 187, 425. 
3 Rees & Ostriker, 1977, MNRAS, 179, 541. 
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Smaller stellar systems: 
Stellar clusters 


While this book focuses on galaxies, there exist smaller gravitating sys- 
tems in which very interesting dynamical mechanisms that are not found 
in galaxies play an important role. Stellar clusters, which inhabit galax- 
ies themselves, come in two flavours: open clusters and globular clusters. 
This chapter explores the properties of these systems and the differ- 
ences with respect to the dynamics of a galaxy. Strong tidal forces are 
especially important, leading to the destruction of lower-mass clusters. 
In addition, the high densities found in the cores of globular clusters 
lead us to explore in more detail the connection between collisionless 
dynamics and thermodynamics. As the velocities of stars in a globular 
cluster evolve towards a Maxwell-Boltzmann distribution, a small but 
sizeable fraction of stars acquires high speeds and leaves the system. 
This process leads to a slow but gradual evaporation. The King mod- 
els are a simple way of describing these systems. In the limit of very 
high densities, a collisional, gravitating structure would collapse into a 
singularity (gravothermal catastrophe). However, we do not find massive 
black holes at the centres of globular clusters. Other mechanisms coun- 
teract this collapse, most notably the formation of tightly bound binaries 
(Heggie’s law). 


8.1 Open and globular clusters 


Stellar clusters are sites of past star formation activity, where gas clouds 
cooled and fragmented. Their dynamical evolution is strongly affected 
by their motion in the Galaxy. The different morphological appear- 
ances of open and globular clusters (figure 8.1) reflect a substantially 
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Figure 8.1 Image of an open cluster (left, Pleiades) and a globu- 
lar cluster (right, M13). DSS-2 images, and inset of M13 taken by 
the Hubble Space Telescope/ACS. (Source: courtesy Space Telescope 
Science Institute.) 


disparate formation scenario. Globular clusters (GCs) are very dense stel- 
lar systems, where up to a million stars can be found within a spherical 
volume of a few parsec in radius. Their stellar populations are overall very 
old and metal poor. The orbits of GCs are not confined to the plane of 
the Galaxy, and are instead found over a large spheroidal volume. These 
properties suggest an early formation process during the collapse of the 
gas component towards the centre of the dark matter halo. Our Milky Way 
galaxy contains about 200 GCs, with half of them located within 10 kpc 
of the Galactic centre, but with a number of them being very distant 
(250 kpc). 


Exercise 8.1 


Consider the motion of a globular cluster of mass m moving along 
a circular orbit around a galaxy with an isothermal potential p(r) = 
v2 /4Gr2, where ve is the (constant) orbital velocity. Dynamical 
friction (see section 3.10) will slow down the orbit. Noting that this 
force acts tangentially along the orbit, find the rate of decrease of 
angular momentum and the time it would take to sink to the centre 
of the galaxy from some initial distance rọ. Assume constant orbital 
velocity throughout. 
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In contrast, open clusters (OCs) comprise younger and more metal 
rich populations, and their orbits are mostly confined to the galactic 
plane. Moreover, they are less dense systems than GCs, and there is an 
age-compactness relation, so that older OCs are more dispersed. These 
properties show that OCs are typical sites of star formation in the Galaxy. 
In the past, about 4.5 Gyr ago, the Sun formed as part of an open cluster. 
Young systems such as the Orion nebula provide a glimpse into the very 
formation of an OC. 


8.2 Internal evolutionary effects 


We will follow a simple argument based on the gravitational energy con- 
tent of stellar clusters to explore the evolution of these systems under the 
effect of internal mechanisms. Starting from the virial theorem, 


GM? 
W=2E = -y —, (8.1) 
R 
and assuming y does not evolve, i.e., the stars in the cluster keep the same 


mass distribution, we find that 
(8.2) 


Let us examine the effects of two important mechanisms. 


Mass loss by stellar winds 


Stellar winds are processes linked to stellar evolution, where the outer 
layers of stars are being pushed away, for instance in massive Wolf-Rayet 
stars, or in the later evolutionary phases (e.g., Asymptotic Giant Branch) 
of lower mass stars. We could also include here the mass ejection from 
supernovae explosions, as massive stars reach their endpoints. If a mass 
dM escapes from the cluster, there is a net change in kinetic energy: 


dT= lam ((v2) — 0) = = ?)dM = —dE. (8.3) 


The last steps involve the virial theorem. Using (v?) = —2E/M, 


dE „dM dR dM 
7M? R zr © RX 1M. (8.4) 
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Consequently, mass loss leads to cluster expansion, by a factor of two or 
more during the lifetime of a cluster. If the mass loss takes place very 
quickly — for example, as a result of young hot stars blowing away the pro- 
tostellar material — then we have an energy Eo = To + Wo = Mo (v2) /2 — 
yGM2 /R, before the event, and E = BT + BW, after the event, where 8 
is the mass fraction remaining in the cluster (we assume the mass loss is 
so rapid that both (v2) and R do not have time to adjust). Hence 


E= p26 —1)Eo. (8.5) 


If half of the mass were lost, we total energy E=0, which means the 
cluster becomes unbound. This is one of the mechanisms that explain why 
so few stars (<1 per cent in the Galactic disc) live in open clusters. 


Core collapse: The gravothermal catastrophe 


Consider a cluster where evolution is self-similar (i.e., the shape and mass 
distribution do not change). Invoking once more the virial theorem, we 
can write: 


E=——~—_, (8.6) 


where y is a dimensionless parameter that depends on the density profile. 
If mass loss is caused only by evaporation, we can neglect the change in 
energy. Hence 


2 2 2 
y GM5 y GM M 
= R = R Da . 8. 
2 Ro 3R > 0 Mo (8.7) 
The mass loss rate can be written 

dM M 
— =-A—, (8.8) 
dt trel 


where A is the mass fraction loss due to evaporation per unit relaxation 
time (i.e., A < 0.01). Here we need to know only the scaling of the relax- 
ation time with respect to the mass and the radius of the cluster. From 
equation 3.23, 


trei x (v?)?/2/ p x v/MR3, 


and using equation 8.7, 


M 7/2 
trel = trel,0 Mo > 
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we can integrate equation 8.8, getting 


2/7 
4 =) , (8.9) 


M(t) =Mo € me 
2 trel,0 


Therefore, the cluster will evaporate in a finite time given by 


7 Atev 2trel 0 
~1lSte~ =~ ~ A0trel.o> 
2 fete ev 7h rel,O 


where the last estimate uses à ~ 0.7%. Turning now to the evolution of 
the stellar density in the cluster, 


The density increases very quickly with time, leading to core collapse. 
Tidal forces will accelerate the evaporation rate. This runaway process 
is inherently caused by the properties of the gravitational force. If we 
model the stars in a cluster as a gas of particles, in analogy with statist- 
ical mechanics, we find that this gas features a negative heat capacity (as 
‘temperature’ increases, the internal energy gets more negative). There- 
fore, as the system ‘loses’ thermal energy, it gets hotter, i.e., the velocity 
dispersion rises. The system will evolve away from an isothermal distri- 
bution, towards a state of higher entropy consisting of an isothermal core 
and a nonisothermal envelope. The core collapses, becoming more tightly 
bound (a cuspy structure), while the halo expands. A singularity may be 
found at the centre. Roughly 20 per cent of globular clusters have central 
density profiles indicative of core collapse. Unsurprisingly, these are the 
clusters with the shortest relaxation times. 


Additional long-term internal evolutionary effects 


Mass segregation: When a system has a shorter relaxation time, stars 
will interact very frequently, exchanging energy towards equipartition. So 
far, we assume all stars have the same mass, leading also to an equiparti- 
tion of the velocity distribution. However, a more realistic scenario would 
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feature a range of stellar masses. In this case, massive stars will transfer 
energy, on average, to lower mass stars, leading to mass segregation, 
where the lighter stars speed up, and the massive ones slow down, sinking 
to the cluster centre. 


Post core-collapse evolution: As density rises during core collapse, the 
number of star-star interactions increases significantly, even leading to 
physical collisions — usually deemed unlikely in the densities expected of 
galaxies. There are typically hundreds of collisions during the lifetime of a 
dense cluster. Blue stragglers are stars that appear on the Main Sequence 
above the turn-off point; i.e., they are more massive than expected from 
the age of the cluster. These stars were created by collisions. Such a high 
density environment can also lead to the formation of a singularity (i.e., 
black hole), formed at the centre of the cluster and expected to be of in- 
termediate mass (~10°Mo). A few possible cases are known so far (e.g., 
M15). The paucity of collapsed cores is likely the result of a mechanism 
that opposes this process: heating by binary systems, an effect we will see 
next. 


Stellar binaries and Heggie’s law: So far, we have considered only indi- 
vidual stars, but they can exist as binary systems, and their internal energy 
can help ‘absorb’ the increasing binding energy as the core collapses. Bin- 
aries can be formed in a primordial way. Pre-collapsed cores contain stars 
formed already as binaries. This is still a relatively unknown quantity, but 
it could be very significant. They could also form via three-body encoun- 
ters, or even in a two-body encounter if tidal interactions are involved, so 
that variations in the internal energy can compensate for the evolution 
from an unbound (positive energy) to a bound system (negative energy). 
Encounters with binaries will change the orbits of stars, acting as a sink 
or a source of gravitational energy. 

Binaries can be defined as soft or hard depending on whether their 
binding energy is less or greater than the mean kinetic energy of indi- 
vidual stars, respectively. Dynamical studies show that as a result of en- 
counters, soft binaries (Ep « mo“) get softer (decreasing |Ep|), and hard 
binaries (Ep >> mo?) get harder (increasing |E,|). This is known as ‘Heg- 
gie’s law’.! This behaviour can be explained via equipartition. Consider a 
soft binary. If the kinetic energy of field stars is greater than the kinetic en- 
ergy of stars in the binary, the internal energy of the binary will be raised, 
becoming softer. This process will eventually lead to the break up of the 
binary, operating over timescales shorter than the relaxation timescale. 
We can envisage an equilibrium state when the rate of binary loss equates 
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the rate of binary formation via three-body encounters. Hence, this pro- 
cess is not important in the evolution of the core. In the opposite case, a 
hard binary has a higher kinetic energy than a typical field star. The inter- 
action will transfer energy to the field star, making the binary more bound 
(harder). This can be an important energy source in dense star clusters, 
and could stop core collapse, and even drive the expansion of the core, at 
the expense of hardening the binary population. The exchange of energy 
with binaries makes X-ray binaries (systems that emit in X-ray because of 
the infall of gas from a main sequence star to a compact object in a close 
binary) more common in the central regions of globular clusters. 


Exercise 8.2 


Let us describe a GC as an isothermal sphere, truncated at r9=5 pc, 
with total mass M=10°Mo. If a mass fraction, f, is in the form of 
close binaries, consisting of two 1 Mo stars orbiting with a sep- 
aration a=10 Ro, determine the value of f for which the binding 
energy in binaries equals the total binding energy of the GC. 


8.3 External effects: Tidal disruption 


Stellar systems are affected by their motion within the host galaxy. Being 
extended systems, they are especially prone to tidal interactions. Tidal 
forces arise from the fact that an extended distribution of mass feels dif- 
ferent gravitational forces in different regions (such is the case with the 
Moon-Earth system). Therefore, stars in the outer regions of a cluster 
are substantially affected either by the potential of the galaxy or by the 
passage of nearby masses. Moreover, for a GC with an eccentric orbit 
outside of the Galactic plane, the passage through the disc produces a 
tidal pulse. We can therefore classify tidal effects as either continuous 
processes caused by the motion within the potential of the galaxy, or 
tidal shocks, where the effect is experimented during a relatively short 
timescale. We will consider these two mechanisms below. 


Steady tidal forces: Jacobi radius 


Consider a cluster orbiting its host galaxy. Tidal forces exerted by the grav- 
itational potential of the galaxy will be able to strip stars from the cluster. 
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Figure 8.2 Simplified version of the steady tidal interaction as a 
globular cluster (GC) orbits around the Galaxy (Glx). 


Therefore, at a given radius (the tidal radius) the number density of stars 
will drop abruptly as their orbits are no longer bound to the cluster. 

Let us explore a simple case, where the cluster (mass m) orbits 
the galaxy (mass M), at a distance D, on a circular orbit, beyond the 
outer edge of the host system. Both galaxy and cluster are spherical. 
The angular speed of the orbital motion (around the common centre of 


mass) is 
o= fo. (8.10) 


The total energy per unit mass of the cluster can be written in terms of the 
effective potential: 

15 is J2 

E= =v + ber) = v + O() +S. 8.11) 

2 eff ( ) 2 ( ) 2r2 ( 
This effective potential — which includes the contribution of angular mo- 
mentum to the energy — introduces a zero-velocity surface, such that 
orbits will never enter the region where € < ®eg(r). Figure 8.2 shows con- 
tours of the effective potential for the simplest case of two orbiting point 
masses. The real thing will be roughly similar. At some critical value of the 
distance, we have the last zero-velocity contour for the cluster. Hence, we 
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should identify the tidal radius at this point (r = rz), such that 


® 
(e) =0. (8.12) 
or Jr 


In our simple model with two point masses we have 


GM Gm a D ) ? Gan 
D-r)? r D3 l+m/M J7 : 


Given that m <M, we have ry «D, and we can expand (D — ry) asa 
Taylor series, and truncate to the lowest order: 


yap(2 y^. (8.14) 


This is the Jacobi limit of the mass m, which provides an estimate of the 
tidal radius of a cluster. There are several aspects that can make ry differ 
from the tidal radius: 


+ The zero velocity surface is not spherical (a single radius is not 
valid). 

+ Noncircular (stellar) orbits can be bound to the cluster, even though 
they probe out tor > ry. 

+ The cluster does not move along a circular orbit. The difference will 
be especially large with highly elongated orbits. 

+ The cluster orbits within the host. 


In any case, the Jacobi limit from equation 8.14 gives values around tens 
of parsec, similar to the cutoff radii found in the surface brightness profiles 
of globular clusters. 


Tidal shocks 


In these encounters — for example, between open clusters and giant mo- 
lecular clouds — the stars are given a sudden jolt. On average, the energies 
will increase. Let us consider the effect on a static star in the cluster 
located close to the centre, after the passage of a cloud (see figure 8.3). 
Av is perpendicular to the direction of motion of the cloud. (This is 
similar to the relaxation time approximation we did in section 3.3.) Sum- 
ming up the effect in all stars in the cluster gives the total change in 
energy: 
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Figure 8.3 Illustration of a tidal shock as a cloud passes near a stellar 
cluster. 


AE= Í Me(Ay)?, (8.15) 


We follow by integrating over all encounters with clouds at all impact 
parameters, to get 


(8.16) 


We can use the virial theorem 2T + W = 0 > E = W/2 and take a generic 
expression of the potential energy as 


GM? 
W=- a 8.1 
VR (8.17) 
where y is, once more, a factor depending on the geometry of the cluster. 
Hence 4 
dE y GM< dR, 
— => —. 8.18 
dt 2 R2 dt Hea 
Comparing equations 8.16 and 8.18, we get 
dR. 8G Mnep RÅ (8.19) 


dt  3yv Piieb M, R2 


neb 


We can relate this equation to the evolution of the density in the cluster: 
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; (8.20) 


which means the cluster density decreases at a constant rate. The life- 
time of the cluster against tidal shocks can be estimated by assuming the 
change in density is equivalent to the density itself (Ap: ~ pc), leading to 


2 
YY Pe Reb 
tshock = 25 


6G Pneb Mneb i 


(8.21) 


proportional to pc. Hence, this effect is very important in open clusters 
(lower density and located on the Galactic plane). Higher density clusters 
will be less likely to be affected by tidal shocks (rather, affected by cluster 
evaporation). 

Globular clusters would be less affected by tidal shocks because of 
their higher densities and orbits away from the plane of the disc. However, 
they can be affected by the passages through the disc, and it is believed 
that low-mass globular clusters were indeed dissolved during the early 
phases of the formation of the Galaxy. 


8.4 Cluster evaporation: King models 


One of the differences between a thermodynamic system and a collision- 
less set of particles moving under gravity is the issue of escape velocity: 
In a relaxed cluster, achieving a Maxwellian distribution in the velocities 
does not imply settling into equilibrium. Stars in the high-velocity tails of 
the distribution will escape from the gravitational potential. A new equi- 
librium, after several relaxation timescales, would now imply a slightly 
lower escape velocity, continuing the evaporation process. The simplest 
approximation to solve this problem is to follow an iterative method, al- 
lowing the cluster to relax to a Maxwellian distribution, then ‘chopping 
off’ the high velocity tails at |v| > Vesce, etc. This introduces the concept 
of a loss rate per relaxation period. The distribution function in velocity 
space is 


1 moe) (v2) 
= /2 w = 
M r)a De BCC act 


From the virial theorem we obtain (v2,.) = 4(v) = 120. So, the fraction 
of stars moving faster than the escape velocity is 
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Figure 8.4 Lowered Gaussian distribution to reconcile a Maxwell- 
Boltzmann profile with the expected ejection of unbound stars from a 
cluster (hatched region). 


F ee is e120? any2dy al > e-* x2dx= 0.007 
= JU = ——— = . . 
(27)3/203 Jy... Vit Js 


A more accurate estimate using the Fokker-Planck equation, assuming 
all stars have the same mass, gives F = 0.008. So, roughly 1 per cent of 
the stars escape from the globular cluster during each relaxation time. 
Therefore we can approximate the distribution function as a lowered Max- 
wellian velocity distribution (see figure 8.4), where the value of the phase 
space distribution function is identically zero for speeds greater than vesc. 
There are two aspects that simplify the dynamical modelling of globular 
clusters: 


+  trelax > tcrossing => the collisionless Boltzmann equation (CBE) can 
be used. 
+  trelax < age of the cluster > clusters are relaxed. 


We need to obtain a self-consistent density-potential pair (0, ®) compat- 
ible with the lowered Gaussian distribution function. For a spherically 
symmetric system there is only one spatial variable, r; the velocity can be 
radial or transverse, (v,,v;); and there are only two integrals of motion: 
Energy (E€ =y7/2+ ) and angular momentum (J =rv;). Therefore, the 
density profile can be written: 


p= f fe. ves y\davedndy, = 2m | PE, Dvdvdv, (8.22) 
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where the second equality derives from Jeans theorem. The other 
equation needed to relate the density and potential is Poisson’s equation 
(equation 2.17), in spherical coordinates: 


1d (e22 


ae dr )=4r6p0, (8.23) 


King models provide a practical representation of the density and poten- 
tial of a globular cluster. They make the following assumptions: 


1. The velocity distribution is isotropic, so that F = F(E). 
Because of external tidal forces, the cluster has a finite radius (i.e., 
the Jacobi radius), where we define the velocity to be zero. 

3. The velocity distribution at r=0 is given by the lowered Gaussian 
approximation, namely: 


—y2 2 —y2 2 
k(e v*/20* _ e Vésc/ 20 ) V <Vese, 


f(r=0,v) = (8.24) 


0 V > Vesec- 


Notice ø is the root mean square dispersion along any of the components 
of the velocity. By defining f (0, v), we can calculate f(r, v) at all radii by 


using ve = —2(r). From this we can derive the density distribution: 


p(r= 25/27 ko 3e%0/7* 
JZ _o2.( [®\ T 27 oy” 
[re /0” erf 52 523 z2 , (8.25) 


and the gravitational potential: 


or)-0=6 | eo as=ane|—= | pods + | posas]. 
o s r Jo 0 


(8.26) 


The integral is solved by parts (f udv=uv— f vdu), choosing u = M(<s) 
and applying the Leibnitz integral rule. These equations are solved by 
choosing a central value for the potential (®9) and integrating the above 
equations outwards until &(R) = 0. We need to solve these equations via 
numerical integration, yielding a family of models characterized by the 
concentration index: 


c=log (+) (8.27) 
0 
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where r;=R is the tidal radius (i.e., the ‘total’ radius of the globular 
cluster), and ro is the King radius, defined by 


902 
= ; 8.28 
ro=ų/ ETET (8.28) 


which corresponds to the radial distance where the projected density of 
the isothermal sphere is 0.5013 of its central value; hence ro is close 
to the definition of the core radius. Typically, the concentration lies in 
the region 0.5 <c < 2.3. King models give acceptable fits to globular and 
open clusters. They can also be applied to model the surface brightness 
profile of dwarf spheroidal and elliptical galaxies. With a knowledge of 
ro and ø? (observed), we can infer the central density, hence obtaining 
the central mass-to-light ratio. This technique is known as ‘core fitting’. 
Most clusters are therefore modelled by three parameters: re, rt, po. In the 
central region, the projected two-dimensional surface mass density (see 
section 2.6) can be approximated by 


Xo 
2: 
1+(£) 


which can be extended to a surface brightness profile if we make an 
assumption about Y = M/L (usually kept constant). This function gives 
a good aproximation to the observed profiles of elliptical galaxies. The 
three-dimensional equivalent is the King profile presented in chapter 2. 
For a galaxy (or an unrelaxed cluster), the assumption of isotropy 


xu(r) = (8.29) 


may not be well justified, and we have to consider anisotropic velocity 
distributions. For instance, the orbits in the outer parts of the cluster may 
be mainly radial. 


Note 


1 Heggie, 1975, MNRAS, 173, 729. 
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9 
Larger stellar systems: Galaxy clusters 


Beyond galaxy scales, an extended structure pervades the Universe in 
the form of filaments over very large distances (210s of Mpc). At the 
confluence of these filaments lie clusters of galaxies. We will briefly ex- 
plore their properties in this final chapter, in relation to all the material 
seen so far. A simple model including hydrostatic equilibrium allows us to 
measure the masses of galaxy clusters, leading to yet another robust proof 
of the presence of dark matter. In addition, gravitational lensing based on 
the deflection of light by the presence of large masses is a useful tech- 
nique. We will present the lensing equation and its application to simple 
mass distributions, showing how one can measure cluster (and galaxy) 
masses. The use of galaxy clusters as cosmological probes is also explored 
here, along with the imprint of large-scale structure on galaxy formation. 


9.1 The most massive structures 


Galaxy clusters are gravitationally bound structures comprising several 
hundred galaxies extended over a radius >1 Mpc. Measurements of 
galaxy redshifts towards a cluster reveal an overdensity at the cluster 
redshift, with a velocity dispersion around ~700-1,000 km s71, Clusters 
shine brightly at X-ray wavelengths (figure 9.1). The X-ray light is 
produced by a large amount of diffuse gas at very high temperatures 
(kT~keV). The morphology of the galaxy distribution and the X-ray emis- 
sion is varied, sometimes featuring substructure and irregular shapes, 
but often showing the smooth morphology expected of a system in equi- 
librium. Clusters are the largest virialized structures in the Universe. In 
the bottom-up hierarchical formation scenario, massive clusters assemble 
at later times (see section 7.3). Therefore, in contrast with galaxies, 
their properties are closer to the expectations from linear theory, making 
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Figure 9.1 Image of galaxy cluster Abell 383 in the optical (left, 
Hubble Space Telescope) and X-ray spectral window (right, Chandra 
Space Telescope). Optical photons originate mostly from the galaxies 
within the cluster, including a few background galaxies distorted in the 
shape of arcs by gravitational lensing (see section 9.3). The X-ray image 
shows the hot diffuse intracluster medium (see section 9.2). 

(Sources: Optical courtesy: NASA/STScI; X-ray courtesy: NASA/CXC/ 
Cinestav/T.Bernal et al.) 


cluster samples valuable cosmological probes. The virial theorem (see 
section 3.9) allows us to determine the typical mass range of a galaxy 
cluster. For an isothermal density profile, we find that 


302R o 2/ R 
M~ —— x6. x 10M, ( ) ( | (9.1) 
G ©\4,000kms-!/ \1Mpe 


Making the simple assumption that a rich cluster includes the equival- 
ent of 100 galaxies with an average luminosity similar to the Milky Way 
(4 x 1010Lo) gives a mass-to-light ratio Y ~ 200Yo, significantly higher 
than any stellar population would produce. Hence the contribution from 
stars to the total mass budget in clusters is negligible. In the next sec- 
tion we will see that X-ray observations allow us to constrain the gas 
mass as well, which, although overwhelmingly higher than the stellar 
mass content, is still short of the masses needed to explain the high velo- 
city dispersion. Clusters are strongly dark matter-dominated systems (see 
figure 1.6), and produce very interesting observational constraints about 
the interplay between dark matter and baryons. 


9.2 X-ray measurements of the cluster mass 
An important tracer of galaxy clusters is the gaseous component 


pervading the cluster in a hot diffuse form. It includes the vast majority 
of the baryon matter content, up to 80-90 per cent by mass. The presence 
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of this so-called intra-cluster medium (ICM) can be explained follow- 
ing the cooling argument invoked in section 7.6. The low density of 
the ICM results in very long cooling timescales, preventing the col- 
lapse and eventual transformation of this gas into stars. It can be de- 
scribed by a relatively simple model, following the equation of state of an 
ideal gas: 

p=nkT = (2) kT, (9.2) 

ump 

where ump is the mean weight of a gas particle, and mp is the proton mass. 
The gas is assumed to be in hydrostatic equilibrium within the gravita- 
tional potential of the underlying dark matter halo. This model allows 
us to derive the mass of a galaxy cluster by observations in the X-ray 
spectral window. We obtain 


dp _ kT dpg i pgk dT GM(<r) 
dr ump dr ` umpdr r2 


Pg, (9.3) 


leading to the following expression: 


M(<r)= (9.4) 


kTr /dlnT l dln pg 
Gump \ dlnr ` dhr J` 


How hot is this gas? By use of the virial theorem, and relating the kinetic 
energy per particle to the thermal energy, we get 


ix = GM 


3kT GM a 
(Fe ) (9.5) 
ump R 


The density of the cluster can be obtained from the spherical col- 
lapse model, adopting the Planck 2015 cosmological parameters (see 
table 7.1), namely: 


kT =0.39 keV M4 no.s (1 +zvR), (9.6) 


where M34 is the cluster mass in units of 10!4Mo, and no.5 = u/0.5. 
Note the virialization of a cluster takes place at relatively low redshift 
zvır ~ 1. Therefore, typical cluster masses yield ICM temperatures in the 
range ~1 keV, i.e., in the X-ray spectral window. The dominant emission 
of the hot ICM gas is produced by thermal bremsstrahlung,! so that the 
emissivity at frequency v is 


J(v) x peT Pe nie, (9.7) 
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A spatially resolved observation of the surface brightness profile in 
X-ray (xJ), including its spectral shape (x dJ/dv), allows us to determ- 
ine the density and temperature profile of the gas, which, via equa- 
tion 9.4, leads to the total mass profile M(<r). A density distribution 
commonly used to model the dark matter in clusters is the King profile 
(see equation 2.34). Noting the gas is in hydrostatic equilibrium (see exer- 
cise 3.6), we understand that the gas density is described by the so-called 


b-model: 
„N2 -368/2 
Pg(T) = po Í + (=) l , (9.8) 


where £ = umpo?/kT is the ratio between the thermal energy of the gas 
and the typical energy of a dark matter particle. By fitting this model to 
the observations, we can derive the gas and total mass profile. A com- 
parison of the total mass with the luminosity in the optical window, i.e., 
contributed by galaxies, gives values of the mass to light ratio of Y > 
200Yo, compatible with our previous virial mass estimate (equation 9.1). 
Moreover, integrating the gas density profile (equation 9.8) with the con- 
straints from X-ray emission produces a mass content that amounts to over 
80 per cent of the total baryon budget in clusters. Note that even in this 
case, the total baryon mass to light ratio would be Yp ~ 50Yo, still lower 
that the expectation from the virial theorem, and unmistakably requiring 
a large fraction of the cluster mass in the form of dark matter. 


9.3 Gravitational lensing 


The gravitational field of a mass distribution distorts spacetime, bending 
the path of photons. This effect is especially strong when the photons 
originate from a distant background source along a very similar line of 
sight to the observer. We can derive the lensing equation by adopting the 
thin lens approximation, a simple, but realistic assumption that the ex- 
tent of the perturber (i.e., the lens) is much smaller than the distance 
between the lens and the source (Dzs) or between the observer and the 
lens (Doz). 

The layout of the lensing problem is shown in figure 9.2. The ob- 
server (O) looks at a distant galaxy (S), located very close to a foreground 
mass distribution that acts as a lens (L). Our measurements are based on 
angular separations on the sky. The angle @ is the observed position of S, 
at point J on the source plane, whereas the true position — which would 
have been the observed position, had the lens not been there — is 6. The 


FUNDAMENTALS OF GALAXY DYNAMICS, FORMATION AND EVOLUTION 


O 
BI of =O 
Source | Lens 
plane plane 
D, S | Da 
D, 1S 


Figure 9.2 Gravitational lensing effect. Light from a background 
source (S) is deflected by a mass (L) between the source and the ob- 
server (O). The deflection angle a causes the observer to image the 
source at position I. The different parameters shown here help us 
derive the lensing equation (9.10). 


deflection angle caused by the effect of the lens is a. For a point mass 
acting as lens, the deflection angle of a light ray hitting the lens plane at 


impact parameter b is 
4GM 
a= ap (9.9) 
We assume that the angles considered here are very small (0, œ, 8 <1), 
so that B =y/Dos and 0 =x/Dos =b/do, and x —y=aDrzs. The lensing 
equation (relating observed and true position of the source) is therefore 


Dıs 4GM Dis 1_ 6% 
Dos c2 DorDos@ 6’ 


4GM Drs 
0 = .11 
í Vy c? DoDos a 


is the Einstein radius of the lens. This is a second order equation in the 
observed position, with solutions 


where 
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B+ ,/B2 +402 
62 — po -02 =0> 0. = 5 E (9.12) 


If 8 =0 - i.e., the background galaxy is exactly positioned along the line 
of sight to the lensing mass — the source will appear as a ring with radius 
Og (Einstein ring). Otherwise, the equation produces two solutions: one 
image ‘inside the ring’ (0 < 0g) and another one ‘outside’. For an extended 
mass distribution with spherical symmetry, the solution is analogous, re- 
placing the mass in the deflection angle (equation 9.9) by the cumulative 
mass contained within the impact parameter, M(< b). 

Figure 9.3 shows an example of the use of the image positions to 
solve for the mass contained inside. Note that some asssumption has 
to be made regarding the mass density profile of the lens. For a case 
close to total alignment (6 < 1), the mass contained within the Einstein 
radius should be independent of this assumption. In general, lensing mass 
estimates tend to have uncertainties that sharply decrease at the equi- 
valent of the Einstein radius. As the true position of the source departs 
outwards (x = 6/8 <1), one of the solutions is very close to the true po- 
sition, 6 = (1 + 2x), and the other one is close to the origin (0@_ = —2Bx) 
and is therefore not visible as it is located at the position of the lens. This 
behaviour defines the weak lensing regime, whereas £ ~ 6g corresponds 
to strong lensing. If we relax the assumption of spherical symmetry for 
the lens mass distribution, more complex solutions are produced.” To put 
in context the expected behaviour in galaxies and clusters, the deflection 
angle in typical units of these systems are 


=i 
a = 0.4 arcsec M b 
1011Mo / \ 10 kpe 


= 40 arcsec ( M ( i 7 (9.13) 
7 10145Mp/\1Mpc/ ’ i 


making strong lensing events available over both galaxy and cluster scales 
— noting that typical ground-based images without adaptive optics can 
reach a spatial resolution ~0.5-1.0 arcsec, and the Hubble Space Tele- 
scope features a typical resolution ~0.1 arcsec. An additional quantity 
often used in gravitational lensing is the critical surface mass density: 


c*Dos 


x, = ——_—_. 
4r GDrsDoL 


(9.14) 


This expression depends only on the geometry of the lens, i.e., the angular 
diametre distance to lens and source and between them, given by the 
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Figure 9.3 Luminosity (dots) and lensing mass profile (shaded area) 
of the Einstein Cross, a lensing system that produces four images of a 
background quasar. Note that the uncertainties in the lensing mass pro- 
duce a typical “butterfly” diagram, where the uncertainty is lowest at 

the position of the Einstein radius, represented by a circle in the inset. 
(Source: Ferreras et al., 2010, MNRAS, 409, L30.) Used by permission of 
Oxford University Press. 


redshifts and the adopted cosmology. £. represents the typical surface 
mass density needed for the presence of a detectable lensing signal, i.e., 
with large enough g. This is especially evident when comparing it with 
the definition of the Einstein radius: 


M 


= ear P 


Xc 


Exercise 9.1 


Consider the strong gravitational lens HE1104-1805%, a galaxy 
at redshift z, = 0.73 that produces two images of a background 
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source at redshift zs = 2.32, with an apparent angular separation 
of 1.1 and 2.1 arcsec away from the galaxy in opposite directions. 
Taking a fiducial ACDM cosmology (h=0.7; Qm,o =0.3), the 
angular diameter distances are: Do, =1498 Mpc; Dos =1690 Mpc; 
and Dzs =910 Mpc. Assuming a point source model for the lens, 
calculate the lensing mass. What would be your answer if you 
adopt an isothermal profile? 


7https://www.cfa.harvard.edu/castles. 


9.4 Clusters and cosmology 


In chapter 7 we traced the evolution of galaxies from the linear regime, 
where the density fluctuations are very small, to the collapse and virializa- 
tion of a structure. We could also appreciate how the standard ACDM (see 
chapter 7) framework leads to a bottom-up formation, so that the most 
massive structures virialize, on average, at later times. Galaxies are sup- 
posed to collapse at rather early cosmic times with respect to the present, 
so that linear theory is not applicable to describe the properties of nearby 
galaxies. In contrast, clusters, or, even better, superclusters, are more 
massive structures whose formation took place at later times, thereby al- 
lowing us to use linear theory to describe their evolution. The distribution 
of clusters and their evolution with redshift thus provides a useful probe of 
the cosmology. For instance, one of the key cosmological parameters that 
describes the amplitude of the power spectrum of density fluctuations is 
og, i.e., the variance (when squared) of the mass fluctuations over a scale 
of 8h—! Mpc, a size that engulfs the typical mass of a galaxy cluster. One of 
the most direct applications of the Press-Schechter methodology (see sec- 
tion 7.4) is to infer this normalization from the number density of clusters 
at present time. Observations of X-ray clusters at low redshift? give a num- 
ber density of M; = 10!5Mo clusters around 3.8 x 1078 Mpe~? within a 
mass interval of 10!+Mo (in this exercise we are assuming h = 0.7 and 
Qm,o = 0.3). We use equation 7.61, noting that the average background 
density at present time is (0) = Perit, 2m,o = 4-1 x 10!°M, Mpc~3. The 
mass variance o(M) can be related to the radius (R) of the spherical 
top-hat filter over which we perform the Press-Schechter smoothing:* 
o (R) =ogR,”, where Rg is the radius in units of 8 h-t Mpc, and y = 0.8. 
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Finally, we get 


1 6 1 x2 
Mpc” 10°**Mo Mpc’ 10**Mo 


n(M;) = 3.8 x 1078 
where x= ôR} /ogV2. The solution to this equation is x~ 2.04. If we 
consider that the linear density contrast for collapse is 6, ~ 1.7 (see equa- 
tion 7.57) and take Rg = 1, we infer a normalization of og = 0.6. This very 
simple exercise gives the right ballpark value (see table 7.1), illustrat- 
ing how clusters can be approximately considered ‘linear structures’ to 
characterize the normalization of the power spectrum.° 


9.5 Environment-related processes 


In addition to the cosmological importance of galaxy clusters as tracers 
of structure growth at late times, they also enable us to probe the vari- 
ous mechanisms of the baryon physics that drives galaxy formation. The 
simple description of an isolated ‘structure unit’ comprising a galaxy living 
in the centre of a virialized dark matter halo, along with an additional dif- 
fuse gas component pervading the halo, breaks down when we take into 
account the complex environment of galaxies over larger scales, with such 
units under constant gravitational interaction within a highly nontrivial 
cosmic web (see figure 1.4). In this context, clusters represent regions 
where galaxies experiment the extreme case of environment-related 
processes. Large surveys in combination with highly detailed numerical 
models of galaxy formation in a cosmological context have allowed us 
to visualize the potential factors that regulate the evolution of galaxies. 
The most important environment-related effects, in no particular 
order, are: 

1. Ram-pressure stripping, suffered by the gas component in galax- 
ies as they fall into a cluster: This process will be more efficient at increas- 
ing infall velocity and at lower surface gas density (i.e., the outer regions 
of galaxies are more prone to stripping). The pressure is exerted by the 
hot intracluster medium. This process is one of the factors that control the 
quenching of star formation in galaxies as they enter clusters, potentially 
transforming spiral galaxies into lenticulars (SO). 

2. Harassment, which is a repetitive process of gravitational interac- 
tion among galaxies at high speed: Like the tidal shock scenario presented 
in section 8.3, each high-speed passage acts as an impulsive force on the 
stars and gas, increasing the energy of the galaxy, whose components 
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Figure 9.4 Variation of the stellar population properties in Early-type 
galaxies according to their local ‘mass’ (parameterized by the velocity 
dispersion, o ) and their environment (split into galaxies in low- and 
high-mass groups). From top to bottom, the age, metallicity, abund- 
ance ratio and dust attenuation are shown with respect to o. At fixed 
galaxy mass, galaxies in more massive halos appear younger, dustier and 
with a lower [a/Fe], representative of a more extended star formation 
history. (Source: La Barbera et al., 2014, MNRAS, 445, 1977.) Used by 
permission of Oxford University Press. 


become less bound. Ina cluster environment, these events will be quite re- 
petitive, leading to the disruption of the galaxies, especially by removing 
the outer, less bound envelopes, and by triggering dynamical instabilities 
that will affect the morphology. 

3. Strangulation: Star formation in a galaxy can be fuelled by the in- 
fall of gas from a diffuse component throughout the halo. As a galaxy falls 
into a cluster, this reservoir may be depleted. Note the difference from ram 
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pressure stripping, which removes the interstellar medium of the galaxy, 
therefore instantaneously truncating the star formation process. The ef- 
fect of strangulation on the star formation rate is more gradual, extending 
over larger timescales. 

4. Mergers and cannibalism: The interaction between two galaxies 
and their host halos may lead to a merger, greatly helped by dynamical 
friction (see section 3.10). Note that a merger process is more effective if 
the relative velocities of the galaxies is comparable to the velocity disper- 
sion of the underlying stellar components. Therefore, the cores of massive 
clusters, with relative velocities ~1,000 km s~! are not fertile ground, but 
lower mass groups and the outskirts of clusters are regions where galaxy 
mergers may be present. An additional merging process, termed ‘canni- 
balism’, can be quite efficient in the central regions of clusters, where 
galaxies will lose energy via dynamical friction, merging with the central, 
most massive galaxy. 

Therefore, if we assume that, to the lowest order, the formation of a 
galaxy is driven by its total mass (i.e., a local property), the combination 
of all these environment-related processes produces an entangled mixture 
of star formation histories that greatly complicate the analysis and extend 
the formation to a global context. In some cases, for instance choosing 
early-type galaxies - whose more homogeneous stellar populations are 
easier to study — it is possible to quantify the role of environment, as shown 
in figure 9.4, where the spectra of early-type galaxies from the Sloan 
Digital Sky Survey were used to constrain the properties of the stellar 
populations. At fixed velocity dispersion (a rough proxy of galaxy mass), 
central galaxies living in more massive halos appear slightly younger, with 
lower [w/Fe], reflecting a more extended star formation history than their 
counterparts in lower density environments. This effect can be explained 
by the contribution of infalling satellite galaxies. 


Notes 


1 Bradt, 2008, Astrophysics processes, Cambridge, section 5.5. 

2 Schneider, Kochanek & Wambsganss, 2006, Gravitational lensing: Strong, weak and 
micro, Springer. 

3 E.g., Reiprich & Bohringer, 2002, ApJ, 567, 716. 

4 Viana & Liddle, 1996, MNRAS, 281, 323. 

5 Allen, Evrard & Mantz, 2011, ARA&A, 49, 409. 
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‘Dr Ferreras strikes gold here with a precisely targeted exposition of the essentials in 
galactic studies. Clearly derived from the extensive experience that he and his colleagues 
at both Oxford and UCL have gathered in teaching this material, it precisely fulfils the 
needs of advanced undergraduate and postgraduate scholars.’ 

— DUNCAN S. MACKAY, Centre for Astrophysics & Planetary Science, University of Kent 


Galaxies, along with their underlying dark matter halos, constitute the building blocks of 
structure in the Universe. Of all fundamental forces, gravity is the dominant one that drives 
the evolution of structures from small density seeds at early times to the galaxies we see 
today. The interactions among myriads of stars, or dark matter particles, in a gravitating 
structure produce a system with fascinating connotations to thermodynamics, with some 
analogies and some fundamental differences. 


Ignacio Ferreras presents a concise introduction to extragalactic astrophysics, with 
emphasis on stellar dynamics, and the growth of density fluctuations in an expanding 
Universe. Additional chapters are devoted to smaller systems (stellar clusters) and larger 
ones (galaxy clusters). Fundamentals of Galaxy Dynamics, Formation and Evolution is 
written for advanced undergraduates and beginning postgraduate students, providing a 
useful tool to get up to speed in a starting research career. Some of the derivations for the 
most important results are presented in detail to enable students appreciate the beauty 

of maths as a tool to understand the workings of galaxies. Each chapter includes a set of 
problems to help the student advance with the material. 


Ignacio Ferreras is a staff astronomer at the Instituto de Astrofisica de Canarias, in 
Tenerife, Spain and holds an honorary professor position at the Physics and Astronomy 
department, UCL. He was academic staff at the Mullard Space Science Laboratory, UCL 
for eleven years. 
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